diff --git a/DESCRIPTION b/DESCRIPTION index e0994fc0..e5271766 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -140,7 +140,6 @@ Collate: 'MultiStudyQtlDataset.R' 'QtlFineMappingResult.R' 'SldscData.R' - 'TupleRangesView.R' 'causalInferencePipeline.R' 'colocPipeline.R' 'qtlSumStats.R' diff --git a/NAMESPACE b/NAMESPACE index 90d26e3f..95e2ca3a 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -4,10 +4,8 @@ S3method(as.data.frame,GwasSumStats) S3method(dplyr::arrange,RangedTupleList) S3method(dplyr::filter,RangedTupleList) S3method(dplyr::group_by,RangedTupleList) -S3method(dplyr::group_by,TupleRangesView) S3method(dplyr::mutate,RangedTupleList) S3method(dplyr::select,RangedTupleList) -S3method(dplyr::select,TupleRangesView) S3method(dplyr::slice,RangedTupleList) S3method(dplyr::summarise,RangedTupleList) S3method(postprocessFinemappingFit,mvsusie) @@ -309,9 +307,12 @@ exportClasses(AnnotationMatrix) exportClasses(ColocBoostResult) exportClasses(ColocResult) exportClasses(ColocResultBase) +exportClasses(CtwasResult) +exportClasses(CtwasResultEntry) exportClasses(FineMappingResultBase) exportClasses(FineMappingRow) exportClasses(GwasFineMappingResult) +exportClasses(GwasSumStats) exportClasses(H2Estimate) exportClasses(LdData) exportClasses(LdEigen) @@ -321,10 +322,10 @@ exportClasses(MashPrior) exportClasses(MultiStudyQtlDataset) exportClasses(QtlDataset) exportClasses(QtlFineMappingResult) +exportClasses(QtlSumStats) exportClasses(RangedTupleList) exportClasses(SldscData) exportClasses(SumStatsBase) -exportClasses(TupleRangesView) exportClasses(TwasWeights) exportClasses(TwasWeightsRow) exportMethods("$") @@ -332,6 +333,7 @@ exportMethods("$<-") exportMethods("[") exportMethods("[[<-") exportMethods(as.data.frame) +exportMethods(bindROWS) exportMethods(colnames) exportMethods(colocboostPipeline) exportMethods(computeLdScores) @@ -492,6 +494,7 @@ importFrom(S4Vectors, "mcols<-", DataFrame, SimpleList, + bindROWS, endoapply, mcols, queryHits, diff --git a/R/AllClasses.R b/R/AllClasses.R index eee8b165..6b553525 100644 --- a/R/AllClasses.R +++ b/R/AllClasses.R @@ -79,7 +79,8 @@ setClassUnion("LdMixtureWeights", c("numeric", "NULL")) # ----------------------------------------------------------------------------- # Shared parent of the QTL and GWAS summary statistics collections. # Concrete subclasses (QtlSumStats, GwasSumStats) inherit from -# RangedTupleList and share the ldSketch / genome / qcInfo slots. Each element +# RangedTupleList and share the ldSketch / qcInfo slots (the genome build +# lives in seqinfo, not a slot). Each element # is one tuple's per-variant GRanges: x[[i]], formerly x$entry[[i]]. # # getZ / getN / getMaf / nSnps are @@ -92,7 +93,8 @@ setClassUnion("LdMixtureWeights", c("numeric", "NULL")) #' @description Virtual base class for QTL and GWAS summary statistics #' collections. Concrete subclasses (\code{QtlSumStats}, \code{GwasSumStats}) #' inherit from \code{\linkS4class{RangedTupleList}} and share the -#' \code{ldSketch} / \code{genome} / \code{qcInfo} slots. +#' \code{ldSketch} / \code{qcInfo} slots, and the genome build in +#' \code{seqinfo()}. #' #' Each element is the per-variant \code{GRanges} of one tuple, so #' \code{x[[i]]} is that tuple's summary statistics and the identity columns @@ -102,7 +104,6 @@ setClassUnion("LdMixtureWeights", c("numeric", "NULL")) #' or \code{NULL}. Optional: LD-free workflows (e.g. mash, which operates #' across conditions per variant) carry \code{NULL}; pipelines that need LD #' validate its presence when they consume the collection. -#' @slot genome Character, genome build label. #' @slot qcInfo A \code{list} recording which QC steps ran. Empty \code{list()} #' on construction; populated by \code{summaryStatsQc()} with a per-step audit #' record (filter names, drop counts, liftover target, RAISS settings, etc.). @@ -115,7 +116,6 @@ setClass( contains = c("VIRTUAL", "RangedTupleList"), representation( ldSketch = "LdSketchOrNULL", - genome = "character", qcInfo = "list" ) ) @@ -176,12 +176,51 @@ setClass( gr[onWindowChrom & IRanges::overlapsAny(gr, win)] } +# The build recorded in seqinfo must be exactly one non-NA value. A missing +# build is an error rather than a default: every downstream liftover / LD +# join keys on it, and silently guessing hg38 is how a mismatched panel gets +# through. Mixed builds mean the parts were never comparable. +# @noRd +.sumStatsCheckGenome <- function(object) { + # seqinfo records the build per SEQLEVEL, so a collection spanning none + # has nowhere to keep one -- whether it has no elements at all or only + # empty ones (a PIP-screened region emptied by summaryStatsQc). That is + # the one case where a missing build is not a defect: there is nothing + # for it to describe. A subset that empties an existing collection keeps + # its seqinfo (see .rtlRebuild), so this exempts only what was built with + # no ranges in the first place. + if (length(GenomeInfoDb::seqlevels(object)) == 0L) { + return(NULL) + } + build <- unique(GenomeInfoDb::genome(object)) + build <- build[!is.na(build)] + if (length(build) == 1L && str_length(build) > 0L) { + return(NULL) + } + if (length(build) == 0L) { + return("no genome build in seqinfo(); set one with genome(x) <- ...") + } + str_c( + "seqinfo() names more than one genome build (", + str_flatten(build, ", "), + ")" + ) +} + #' @rdname getGenome #' @examples #' data(qtlSumStatsExample) #' getGenome(qtlSumStatsExample) #' @export -setMethod("getGenome", "SumStatsBase", function(x, ...) x@genome) +setMethod("getGenome", "SumStatsBase", function(x, ...) { + # The build lives in seqinfo, exactly as it does on LdStatistic: a + # GRangesList already has somewhere to keep it, and a parallel `genome` + # slot went stale against it (getGenome() said hg19 while genome(x) said + # NA, so every Bioconductor path that reads genome(x) saw nothing). + build <- unique(GenomeInfoDb::genome(x)) + build <- build[!is.na(build)] + if (length(build) == 0L) NA_character_ else build[[1L]] +}) #' @rdname getQcInfo #' @examples diff --git a/R/CtwasResult.R b/R/CtwasResult.R index 3743b3d2..b14e89a0 100644 --- a/R/CtwasResult.R +++ b/R/CtwasResult.R @@ -14,6 +14,20 @@ #' @include AllGenerics.R tupleSelectors.R NULL +#' @title cTWAS Result Collection +#' @description S4 collection of cTWAS runs keyed by the identity tuple +#' \code{(gwasStudy, study, context, method)}. Each row holds a +#' \code{\linkS4class{CtwasResultEntry}} payload -- fine-mapping posteriors, +#' the jointly-estimated group priors, and region metadata -- for one run. +#' @details Unlike the QTL family, \code{trait} is not part of the key: a cTWAS +#' run is multi-gene, so genes live inside the payload. The optional +#' \code{jointStudies} / \code{jointContexts} columns tag rows born from a +#' multi-study or multi-context run and participate in the uniqueness key, +#' exactly as in the \code{\linkS4class{TwasWeights}} and fine-mapping +#' families. +#' @seealso \code{\link{CtwasResult}} for the constructor and +#' \code{\link{ctwasPipeline}} for the pipeline that builds one. +#' @export setClass("CtwasResult", contains = "DFrame", validity = function(object) { .validateCtwasResult(object) }) diff --git a/R/CtwasResultEntry.R b/R/CtwasResultEntry.R index 30638f92..a296af47 100644 --- a/R/CtwasResultEntry.R +++ b/R/CtwasResultEntry.R @@ -10,6 +10,23 @@ #' @include AllGenerics.R NULL +#' @title cTWAS Per-Run Payload +#' @description Per-run cTWAS payload: the fine-mapping posterior table, the +#' full per-effect susie alpha table, the jointly-estimated group prior(s), +#' and per-region metadata. One entry sits in every row of a +#' \code{\linkS4class{CtwasResult}} collection. +#' @slot finemap The per-gene (and, when SNPs are retained, per-SNP) posterior +#' summary table (\code{ctwas::finemap_regions} \code{finemap_res} shape), +#' or \code{NULL}. +#' @slot susieAlpha The per-effect susie alpha table +#' (\code{ctwas::finemap_regions} \code{susie_alpha_res} shape) -- the +#' fuller cTWAS output retained so the raw run is reconstructable, or +#' \code{NULL}. +#' @slot param The estimated \code{group_prior} / \code{group_prior_var} for +#' this run, or \code{NULL}. +#' @slot regionInfo Per-region metadata, or \code{NULL}. +#' @seealso \code{\link{CtwasResultEntry}} for the constructor. +#' @export setClass( "CtwasResultEntry", representation( diff --git a/R/QtlDataset.R b/R/QtlDataset.R index d28f73ce..cda5972a 100644 --- a/R/QtlDataset.R +++ b/R/QtlDataset.R @@ -46,7 +46,6 @@ NULL #' @slot study Character (length 1). Study identifier; used in collection #' classes to tag downstream \code{FineMappingResult} / \code{TwasWeights} #' entries. -#' @slot genotypes The genotype source for lazy access to dosages. #' The \code{genotype} experiment's assay reads through this handle; the #' extraction accessors read it directly, so that QC can be applied per #' block. @@ -78,7 +77,6 @@ setClass( contains = "MultiAssayExperiment", representation( study = "character", - genotypes = "GenotypeHandle", scaleResiduals = "logical", mafCutoff = "numeric", macCutoff = "numeric", @@ -403,7 +401,6 @@ QtlDataset <- function( "QtlDataset", .qtlRestrictSamples(mae, keepSamples), study = as.character(study), - genotypes = handle, scaleResiduals = isTRUE(scaleResiduals), mafCutoff = as.numeric(mafCutoff), macCutoff = as.numeric(macCutoff), @@ -575,19 +572,15 @@ setMethod("longForm", "QtlDataset", function(object, ..., genotype = FALSE) { x[, which(is_in(ids, as.character(keepSamples))), ] } -# Replace the genotype handle, rebuilding the genotype experiment along with -# it. The handle is deliberately held twice -- once as a slot, once inside -# the assay's seed -- because that is what lets the assay read lazily. Moving -# one without the other would leave the dosages describing a different panel -# from the one the extraction accessors read, so nothing may set the slot -# directly. +# Replace the genotype handle by rebuilding the genotype experiment around +# it. The handle lives in exactly one place -- the assay's seed -- so there +# is no second copy to keep in step; getGenotypeHandle() reads it back. # @noRd .qtlWithGenotypeHandle <- function(x, handle) { exps <- MultiAssayExperiment::experiments(x) gCov <- .qtlColDataMatrix(exps[[.QTL_GENO_EXPERIMENT]]) exps[[.QTL_GENO_EXPERIMENT]] <- .genotypeExperiment(handle, gCov) MultiAssayExperiment::experiments(x) <- exps - x@genotypes <- handle validObject(x) x } @@ -673,7 +666,15 @@ setMethod("getScaleResiduals", "QtlDataset", function(x) x@scaleResiduals) #' @rdname getGenotypeHandle #' @keywords internal -setMethod("getGenotypeHandle", "QtlDataset", function(x) x@genotypes) +setMethod("getGenotypeHandle", "QtlDataset", function(x) { + # Derived, not stored. The handle already lives inside the genotype + # assay's seed -- that is what lets the dosages read lazily -- so a + # parallel slot was a second copy that had to be kept in step by hand. + # Reading it back removes the invariant instead of policing it. + .ldSketchHandle( + MultiAssayExperiment::experiments(x)[[.QTL_GENO_EXPERIMENT]] + ) +}) #' @rdname qtlDatasetFilters #' @export @@ -2069,8 +2070,9 @@ setMethod("show", "QtlDataset", function(object) { .trim = FALSE )) cat(glue(" {totalTraits} unique traits across contexts\n", .trim = FALSE)) + gh <- getGenotypeHandle(object) cat(glue( - " Genotypes: {object@genotypes@format} @ {object@genotypes@path}\n", + " Genotypes: {getFormat(gh)} @ {getPath(gh)}\n", .trim = FALSE )) cat(glue( diff --git a/R/RangedTupleList.R b/R/RangedTupleList.R index 69691c3c..f34abc45 100644 --- a/R/RangedTupleList.R +++ b/R/RangedTupleList.R @@ -132,12 +132,43 @@ methods::setValidity("RangedTupleList", function(object) { # drops the subclass's slots -- verified -- which is exactly the desync this # class exists to prevent. # @noRd +# Put `x`'s seqinfo back onto a rebuilt GRangesList. Only the seqlevels the +# rebuild actually spans can be kept, so this restores the genome (and any +# lengths) for those, and supplies the whole seqinfo when the rebuild spans +# nothing. +# @noRd +.rtlRestoreSeqinfo <- function(grl, x) { + si <- GenomeInfoDb::seqinfo(x) + if (length(GenomeInfoDb::seqlevels(si)) == 0L) { + return(grl) + } + if (length(GenomeInfoDb::seqlevels(grl)) == 0L) { + GenomeInfoDb::seqinfo(grl) <- si + return(grl) + } + keepLevels <- intersect( + GenomeInfoDb::seqlevels(si), + GenomeInfoDb::seqlevels(grl) + ) + GenomeInfoDb::seqinfo( + grl, + new2old = match(keepLevels, GenomeInfoDb::seqlevels(grl)) + ) <- si[keepLevels] + grl +} + .rtlRebuild <- function(x, elements, keep) { # mcols and slots are attached BEFORE new(), not after: new() validates # during initialize(), and a subclass's validity method reads its identity # columns and slots. Building the object bare and filling it in afterwards # trips that check on the way past. grl <- GenomicRanges::GRangesList(elements) + # seqinfo is collection-level state, exactly like the slots below: it + # carries the genome build. A rebuild from bare elements starts with the + # seqlevels those elements happen to span -- none at all when everything + # was dropped -- so the original seqinfo is merged back in, or subsetting + # to nothing would silently discard the build. + grl <- .rtlRestoreSeqinfo(grl, x) md <- mcols(x, use.names = FALSE) if (!is.null(md)) { mcols(grl) <- md[keep, , drop = FALSE] @@ -563,3 +594,313 @@ setMethod("subsetRegion", "RangedTupleList", function(x, region, ...) { } g[onWindowChrom & IRanges::overlapsAny(g, win)] } + + +# ============================================================================= +# Combining collections +# ============================================================================= + +#' @rdname RangedTupleList-methods +#' @importFrom S4Vectors bindROWS +#' @export +setMethod( + "bindROWS", + "RangedTupleList", + function( + x, + objects = list(), + use.names = TRUE, + ignore.mcols = FALSE, + check = TRUE + ) { + # bindROWS() is the one hook `c()` and `append()` share, so defining it + # here fixes both. The inherited CompressedList version rbinds the mcols, + # which requires identical columns and so fails whenever the parts came + # from runs with different optional columns (traitPos / jointContexts / + # blockId); it also carries the FIRST part's collection-level slots + # silently, making the result depend on argument order. + .combineTupleCollections(compact(c(list(x), objects)), NULL, "c") + } +) + + +# ============================================================================= +# The plyranges / dplyr bridge for RangedTupleList collections +# +# plyranges has no GRangesList support at all: its verbs dispatch on +# GenomicRanges / Ranges, and a plain CompressedGRangesList fails exactly as +# these collections do. That is not a deficiency in the collections, and it is +# not fixable by changing them. +# +# The bridge is a round trip through a flat GRanges instead: flatten the +# collection (broadcasting its identity tuple onto every range), let plyranges +# work on the shape it already supports, and nest the result back into the +# collection's elements. Nothing about the collections has to change, and no +# wrapper class is needed -- plyranges acts on a real GRanges natively. +# ============================================================================= + +# ----------------------------------------------------------------------------- +# Conversions +# ----------------------------------------------------------------------------- + +#' Flatten a tuple collection to a plain GRanges +#' +#' Concatenates every element's ranges and broadcasts the collection's identity +#' tuple onto each range, so a plyranges predicate can mix tuple and per-range +#' columns: \code{filter(x, .context == "blood" & pip > 0.9)}. +#' +#' Broadcast columns are prefixed with a dot -- \code{.study}, +#' \code{.context}, \code{.trait}, \code{.method} -- following the +#' tidySummarizedExperiment convention for framework-injected columns +#' (\code{.sample}, \code{.feature}). The prefix is not cosmetic: a +#' collection's identity column can collide with a per-range column of the same +#' name carrying different information. On +#' \code{\link{gwasFineMappingExample}} the outer \code{method} is +#' \code{"susie"} while the per-range \code{method} is \code{"susieRss"} -- +#' the collection's label against the fitter actually used. Prefixing keeps +#' both, and keeps the name stable so a predicate does not change meaning +#' between objects. +#' +#' Non-atomic \code{mcols} columns -- the per-element \code{susieFit} / +#' \code{cvResult} payloads -- are dropped. They describe an element, not a +#' range, so there is no row to broadcast them onto. Use the accessors +#' (\code{\link{getSusieFit}}, \code{\link{getCvResult}}) for those. +#' +#' @param x A \code{RangedTupleList}. +#' @return A \code{GRanges}. +#' @examples +#' data(qtlFineMappingExample) +#' flat <- flattenTupleRanges(qtlFineMappingExample) +#' head(names(S4Vectors::mcols(flat))) +#' @export +flattenTupleRanges <- function(x) { + if (!methods::is(x, "RangedTupleList")) { + msg <- glue( + "`x` must be a RangedTupleList (got {class(x)[[1L]]})." + ) + abort(msg) + } + gr <- .rtlGatherElements(x, seq_len(length(x))) + md <- mcols(x, use.names = FALSE) + if (is.null(md) || ncol(md) == 0L || length(gr) == 0L) { + return(gr) + } + tupleCols <- names(md)[map_lgl(as.list(md), is.atomic)] + reps <- rep(seq_len(length(x)), lengths(x)) + for (nm in tupleCols) { + mcols(gr)[[.rtlDotName(nm)]] <- md[[nm]][reps] + } + gr +} + +# The broadcast name for an identity column. Dotted so it cannot collide with +# a per-range column, and so the name is the same whatever the object holds. +# @noRd +.rtlDotName <- function(nm) { + str_c(".", nm) +} + +#' Re-nest a flattened GRanges into a tuple collection +#' +#' The inverse of \code{\link{flattenTupleRanges}}: ranges are grouped by the +#' identity tuple they carry and returned to \code{template}'s elements, so the +#' collection's class, its per-element payloads and its collection-level slots +#' all survive a plyranges round trip. +#' +#' A tuple with no surviving ranges becomes an empty element rather than being +#' dropped, so the collection keeps its shape and its metadata stays aligned. +#' +#' @param flat A \code{GRanges} produced by \code{flattenTupleRanges} (possibly +#' filtered or mutated). +#' @param template The collection it came from, supplying the tuple grid, the +#' payload columns and the slots. +#' @return An object of \code{template}'s class. +#' @examples +#' data(qtlFineMappingExample) +#' flat <- flattenTupleRanges(qtlFineMappingExample) +#' nestTupleRanges(flat, qtlFineMappingExample) +#' @export +nestTupleRanges <- function(flat, template) { + if (!methods::is(template, "RangedTupleList")) { + msg <- glue( + "`template` must be a RangedTupleList (got ", + "{class(template)[[1L]]})." + ) + abort(msg) + } + if (!methods::is(flat, "GRanges")) { + msg <- glue("`flat` must be a GRanges (got {class(flat)[[1L]]}).") + abort(msg) + } + keyCols <- .rtlTupleKeyCols(template) + dotted <- map_chr(keyCols, .rtlDotName) + wanted <- .rtlTupleKeys( + mcols(template, use.names = FALSE), + keyCols, + n = length(template) + ) + have <- .rtlTupleKeys( + mcols(flat, use.names = FALSE), + dotted, + n = length(flat) + ) + elements <- map(wanted, .rtlPickByKey, flat = flat, have = have) + .rtlRebuild(template, elements, seq_len(length(template))) +} + +# The identity columns to group by: the atomic mcols the flattener broadcasts. +# @noRd +.rtlTupleKeyCols <- function(x) { + md <- mcols(x, use.names = FALSE) + if (is.null(md)) { + return(character(0)) + } + names(md)[map_lgl(as.list(md), is.atomic)] +} + +# One key string per row. With no identity columns every range belongs to the +# single element, which is what a one-row collection means. +# @noRd +.rtlTupleKeys <- function(md, keyCols, n) { + present <- intersect(keyCols, colnames(md)) + if (length(present) == 0L) { + return(rep("", n)) + } + exec(str_c, !!!map(present, .rtlKeyPart, md = md), sep = "\r") +} + +# NA is mapped to a sentinel rather than left alone: str_c() propagates NA, so +# a single NA-valued identity column (varY is NA_real_ on a z-score collection) +# would turn every key into NA and the subsequent `have == key` into a logical +# subscript full of NAs. +# @noRd +.rtlKeyPart <- function(nm, md) { + v <- as.character(md[[nm]]) + if_else(is.na(v), "\u0001NA", v) +} + +# @noRd +.rtlPickByKey <- function(key, flat, have) { + flat[have == key] +} + + +# ----------------------------------------------------------------------------- +# Verbs on the collections themselves +# ----------------------------------------------------------------------------- + +# S3 dispatch reaches a method registered on an S4 VIRTUAL base, so one set of +# methods here serves every RangedTupleList subclass -- the fine-mapping, +# sumstats, TWAS-weight and coloc collections alike. +# +# Shape-preserving verbs flatten, delegate and re-nest, so the caller gets the +# collection back. Reducing verbs return what the verb produces, because a +# summary has no per-element shape to nest into. + +# These verbs delegate to plyranges' GRanges methods, which only exist once +# plyranges' namespace is loaded. Without this the failure is an opaque "no +# applicable method for 'filter' applied to an object of class GRanges" -- +# pointing at the flattened form rather than at the missing package. +# requireNamespace() both checks and loads, so the check is also the fix. +# @noRd +.rtlRequirePlyranges <- function(verb) { + if (!requireNamespace("plyranges", quietly = TRUE)) { + msg <- glue( + "`{verb}()` on a tuple collection needs the plyranges package; ", + "install it, or work on flattenTupleRanges(x) directly." + ) + abort(msg) + } + invisible(NULL) +} + +# `...` is forwarded directly rather than captured with enquos() and spliced: +# splicing hands the verb quosure OBJECTS instead of expressions to evaluate, +# and plyranges then fails with "Argument to filter condition must evaluate to +# a logical vector". + +#' @exportS3Method dplyr::filter +filter.RangedTupleList <- function(.data, ...) { + .rtlRequirePlyranges("filter") + out <- dplyr::filter(flattenTupleRanges(.data), ...) + nestTupleRanges(out, .data) +} + +#' @exportS3Method dplyr::mutate +mutate.RangedTupleList <- function(.data, ...) { + .rtlRequirePlyranges("mutate") + out <- dplyr::mutate(flattenTupleRanges(.data), ...) + nestTupleRanges(out, .data) +} + +#' @exportS3Method dplyr::select +select.RangedTupleList <- function(.data, ...) { + .rtlRequirePlyranges("select") + # The identity columns are kept whatever the selection, or the result + # could not be nested back into its elements. + out <- dplyr::select(flattenTupleRanges(.data), ...) + nestTupleRanges( + .rtlRestoreKeys(out, flattenTupleRanges(.data), .data), + .data + ) +} + +# Takes n ranges FROM EACH element, matching arrange()'s within-element rule. +# +# Unlike the other verbs this does not flatten: slicing the flat set would take +# n ranges in TOTAL, which for a many-element collection empties all but the +# first. Applying the slice per element is both the right semantics and simpler +# than reconstructing the grouping after a flatten. +#' @exportS3Method dplyr::slice +slice.RangedTupleList <- function(.data, ...) { + .rtlRequirePlyranges("slice") + elements <- map(as.list(.data), .rtlSliceOne, ...) + .rtlRebuild(.data, elements, seq_len(length(.data))) +} + +# @noRd +.rtlSliceOne <- function(g, ...) { + dplyr::slice(g, ...) +} + +# Orders WITHIN each element. A global sort of the flattened set followed by +# partitioning on the identity tuple leaves each element internally sorted -- +# the two are the same thing for the partitioned result -- so no grouping is +# needed here. Ordering ACROSS elements would permute the elements themselves, +# which the identity tuple cannot express. +#' @exportS3Method dplyr::arrange +arrange.RangedTupleList <- function(.data, ...) { + .rtlRequirePlyranges("arrange") + out <- dplyr::arrange(flattenTupleRanges(.data), ...) + nestTupleRanges(out, .data) +} + +#' @exportS3Method dplyr::summarise +summarise.RangedTupleList <- function(.data, ...) { + .rtlRequirePlyranges("summarise") + dplyr::summarise(flattenTupleRanges(.data), ...) +} + +# No count() method: plyranges defines none for GRanges either, so there is +# nothing to delegate to. Use summarise(group_by(x, ...), n = n()). + +#' @exportS3Method dplyr::group_by +group_by.RangedTupleList <- function(.data, ...) { + .rtlRequirePlyranges("group_by") + dplyr::group_by(flattenTupleRanges(.data), ...) +} + +# Put back any identity column a selection dropped. Without them the ranges +# carry no tuple and every one would fall into the first element. +# @noRd +.rtlRestoreKeys <- function(out, flat, template) { + dotted <- map_chr(.rtlTupleKeyCols(template), .rtlDotName) + missing <- setdiff( + intersect(dotted, colnames(mcols(flat))), + colnames(mcols(out)) + ) + for (nm in missing) { + mcols(out)[[nm]] <- mcols(flat)[[nm]] + } + out +} diff --git a/R/TupleRangesView.R b/R/TupleRangesView.R deleted file mode 100644 index 34f4d65f..00000000 --- a/R/TupleRangesView.R +++ /dev/null @@ -1,383 +0,0 @@ -# ============================================================================= -# TupleRangesView -- the plyranges bridge for RangedTupleList collections -# -# plyranges has no GRangesList support at all: its verbs dispatch on -# GenomicRanges / Ranges, and a plain CompressedGRangesList fails exactly as -# these collections do. That is not a deficiency in the collections, and it is -# not fixable by changing them. -# -# It IS fixable through the extension point GenomicRanges provides: -# `DelegatingGenomicRanges` is a virtual class that wraps a flat `delegate` and -# *is* a GenomicRanges, and plyranges ships methods for it. A subclass supplies -# a `[` hook and inherits the verb surface instead of reimplementing it. -# -# The conversion is the right shape rather than a workaround: plyranges' own -# `GroupedGenomicRanges` (slots group_keys / group_indices / n / delegate) -# represents grouping as flat ranges plus group metadata, which is exactly a -# RangedTupleList transposed. Same information, other representation. -# ============================================================================= - -#' Flat view of a tuple collection for plyranges verbs -#' -#' A \code{\link[GenomicRanges]{GenomicRanges}} view over a -#' \code{RangedTupleList}: every element's ranges concatenated, with the -#' collection's identity tuple broadcast onto each range. Because it is a -#' \code{GenomicRanges}, plyranges verbs act on it directly. -#' -#' Users do not normally build one. The verb methods on \code{RangedTupleList} -#' convert to a view, delegate, and convert back. -#' -#' @slot delegate The flattened \code{GRanges}. -#' @name TupleRangesView-class -#' @keywords internal -#' @export -setClass("TupleRangesView", contains = "DelegatingGenomicRanges") - -# Build a view over `flat`. -# -# DelegatingGenomicRanges carries its OWN elementMetadata alongside the -# delegate's, and validity requires it to be parallel to the object. Supplying -# only `delegate=` leaves it zero-row against a longer object and the class -# rejects itself with "'mcols(x)' is not parallel to 'x'", so it is set here. -# @noRd -.tupleRangesView <- function(flat) { - methods::new( - "TupleRangesView", - delegate = flat, - elementMetadata = S4Vectors::make_zero_col_DFrame(length(flat)) - ) -} - -# The flat GRanges a view wraps. Internal accessor: TupleRangesView is not -# exported, so this adds no public surface and keeps the `@` in the class file. -# @noRd -.trvDelegate <- function(x) x@delegate - -# The one hook a DelegatingGenomicRanges subclass has to supply. Without it, -# plyranges' subsetting verbs (slice, filter_by_overlaps, join_overlap_*) fail -# with "subscript is a NSBS object that is incompatible with the current -# subsetting operation". -#' @rdname TupleRangesView-class -#' @export -setMethod("[", "TupleRangesView", function(x, i, j, ..., drop = TRUE) { - if (!missing(j)) { - abort("two-argument `[` is not supported on a TupleRangesView.") - } - .tupleRangesView(.trvDelegate(x)[i]) -}) - -#' @rdname show-methods -#' @export -setMethod("show", "TupleRangesView", function(object) { - cat(glue( - "TupleRangesView: {length(.trvDelegate(object))} ranges, ", - "{ncol(mcols(.trvDelegate(object)))} metadata column(s)\n", - .trim = FALSE - )) - invisible(NULL) -}) - -# ----------------------------------------------------------------------------- -# Conversions -# ----------------------------------------------------------------------------- - -#' Flatten a tuple collection to a plain GRanges -#' -#' Concatenates every element's ranges and broadcasts the collection's identity -#' tuple onto each range, so a plyranges predicate can mix tuple and per-range -#' columns: \code{filter(x, .context == "blood" & pip > 0.9)}. -#' -#' Broadcast columns are prefixed with a dot -- \code{.study}, -#' \code{.context}, \code{.trait}, \code{.method} -- following the -#' tidySummarizedExperiment convention for framework-injected columns -#' (\code{.sample}, \code{.feature}). The prefix is not cosmetic: a -#' collection's identity column can collide with a per-range column of the same -#' name carrying different information. On -#' \code{\link{gwasFineMappingExample}} the outer \code{method} is -#' \code{"susie"} while the per-range \code{method} is \code{"susieRss"} -- -#' the collection's label against the fitter actually used. Prefixing keeps -#' both, and keeps the name stable so a predicate does not change meaning -#' between objects. -#' -#' Non-atomic \code{mcols} columns -- the per-element \code{susieFit} / -#' \code{cvResult} payloads -- are dropped. They describe an element, not a -#' range, so there is no row to broadcast them onto. Use the accessors -#' (\code{\link{getSusieFit}}, \code{\link{getCvResult}}) for those. -#' -#' @param x A \code{RangedTupleList}. -#' @return A \code{GRanges}. -#' @examples -#' data(qtlFineMappingExample) -#' flat <- flattenTupleRanges(qtlFineMappingExample) -#' head(names(S4Vectors::mcols(flat))) -#' @export -flattenTupleRanges <- function(x) { - if (!methods::is(x, "RangedTupleList")) { - msg <- glue( - "`x` must be a RangedTupleList (got {class(x)[[1L]]})." - ) - abort(msg) - } - gr <- .rtlGatherElements(x, seq_len(length(x))) - md <- mcols(x, use.names = FALSE) - if (is.null(md) || ncol(md) == 0L || length(gr) == 0L) { - return(gr) - } - tupleCols <- names(md)[map_lgl(as.list(md), is.atomic)] - reps <- rep(seq_len(length(x)), lengths(x)) - for (nm in tupleCols) { - mcols(gr)[[.rtlDotName(nm)]] <- md[[nm]][reps] - } - gr -} - -# The broadcast name for an identity column. Dotted so it cannot collide with -# a per-range column, and so the name is the same whatever the object holds. -# @noRd -.rtlDotName <- function(nm) { - str_c(".", nm) -} - -#' Re-nest a flattened GRanges into a tuple collection -#' -#' The inverse of \code{\link{flattenTupleRanges}}: ranges are grouped by the -#' identity tuple they carry and returned to \code{template}'s elements, so the -#' collection's class, its per-element payloads and its collection-level slots -#' all survive a plyranges round trip. -#' -#' A tuple with no surviving ranges becomes an empty element rather than being -#' dropped, so the collection keeps its shape and its metadata stays aligned. -#' -#' @param flat A \code{GRanges} produced by \code{flattenTupleRanges} (possibly -#' filtered or mutated). -#' @param template The collection it came from, supplying the tuple grid, the -#' payload columns and the slots. -#' @return An object of \code{template}'s class. -#' @examples -#' data(qtlFineMappingExample) -#' flat <- flattenTupleRanges(qtlFineMappingExample) -#' nestTupleRanges(flat, qtlFineMappingExample) -#' @export -nestTupleRanges <- function(flat, template) { - if (!methods::is(template, "RangedTupleList")) { - msg <- glue( - "`template` must be a RangedTupleList (got ", - "{class(template)[[1L]]})." - ) - abort(msg) - } - if (!methods::is(flat, "GRanges")) { - msg <- glue("`flat` must be a GRanges (got {class(flat)[[1L]]}).") - abort(msg) - } - keyCols <- .rtlTupleKeyCols(template) - dotted <- map_chr(keyCols, .rtlDotName) - wanted <- .rtlTupleKeys( - mcols(template, use.names = FALSE), - keyCols, - n = length(template) - ) - have <- .rtlTupleKeys( - mcols(flat, use.names = FALSE), - dotted, - n = length(flat) - ) - elements <- map(wanted, .rtlPickByKey, flat = flat, have = have) - .rtlRebuild(template, elements, seq_len(length(template))) -} - -# The identity columns to group by: the atomic mcols the flattener broadcasts. -# @noRd -.rtlTupleKeyCols <- function(x) { - md <- mcols(x, use.names = FALSE) - if (is.null(md)) { - return(character(0)) - } - names(md)[map_lgl(as.list(md), is.atomic)] -} - -# One key string per row. With no identity columns every range belongs to the -# single element, which is what a one-row collection means. -# @noRd -.rtlTupleKeys <- function(md, keyCols, n) { - present <- intersect(keyCols, colnames(md)) - if (length(present) == 0L) { - return(rep("", n)) - } - exec(str_c, !!!map(present, .rtlKeyPart, md = md), sep = "\r") -} - -# NA is mapped to a sentinel rather than left alone: str_c() propagates NA, so -# a single NA-valued identity column (varY is NA_real_ on a z-score collection) -# would turn every key into NA and the subsequent `have == key` into a logical -# subscript full of NAs. -# @noRd -.rtlKeyPart <- function(nm, md) { - v <- as.character(md[[nm]]) - if_else(is.na(v), "\u0001NA", v) -} - -# @noRd -.rtlPickByKey <- function(key, flat, have) { - flat[have == key] -} - - -# ----------------------------------------------------------------------------- -# Patches for the two verbs plyranges' DelegatingGenomicRanges support misses -# ----------------------------------------------------------------------------- - -# Both work on a plain GRanges and fail on any DelegatingGenomicRanges: -# `select` reports "Cannot select/rename the following columns: seqnames, -# start, end, width, strand", and `group_by` -- for which plyranges defines no -# DelegatingGenomicRanges method at all -- reports "Can't select columns with -# `dots`". Both are upstream gaps, not something about this class. -# -# The fix in each case is to run the verb on the delegate, where plyranges -# works, and re-wrap. - -#' @exportS3Method dplyr::select -select.TupleRangesView <- function(.data, ..., .drop_ranges = FALSE) { - out <- dplyr::select( - .trvDelegate(.data), - ..., - .drop_ranges = .drop_ranges - ) - if (isTRUE(.drop_ranges)) { - return(out) - } - .tupleRangesView(out) -} - -# Not re-wrapped: grouping is plyranges' own representation, and the result is -# a GroupedGenomicRanges rather than a view of a collection. -#' @exportS3Method dplyr::group_by -group_by.TupleRangesView <- function(.data, ..., .add = FALSE) { - dplyr::group_by(.trvDelegate(.data), ...) -} - -# ----------------------------------------------------------------------------- -# Verbs on the collections themselves -# ----------------------------------------------------------------------------- - -# S3 dispatch reaches a method registered on an S4 VIRTUAL base, so one set of -# methods here serves every RangedTupleList subclass -- the fine-mapping, -# sumstats, TWAS-weight and coloc collections alike. -# -# Shape-preserving verbs flatten, delegate and re-nest, so the caller gets the -# collection back. Reducing verbs return what the verb produces, because a -# summary has no per-element shape to nest into. - -# These verbs delegate to plyranges' GRanges methods, which only exist once -# plyranges' namespace is loaded. Without this the failure is an opaque "no -# applicable method for 'filter' applied to an object of class GRanges" -- -# pointing at the flattened form rather than at the missing package. -# requireNamespace() both checks and loads, so the check is also the fix. -# @noRd -.rtlRequirePlyranges <- function(verb) { - if (!requireNamespace("plyranges", quietly = TRUE)) { - msg <- glue( - "`{verb}()` on a tuple collection needs the plyranges package; ", - "install it, or work on flattenTupleRanges(x) directly." - ) - abort(msg) - } - invisible(NULL) -} - -# Unwrap a view back to its delegate; pass a plain GRanges through. -# @noRd -.rtlAsRanges <- function(out) { - if (methods::is(out, "TupleRangesView")) .trvDelegate(out) else out -} - -# `...` is forwarded directly rather than captured with enquos() and spliced: -# splicing hands the verb quosure OBJECTS instead of expressions to evaluate, -# and plyranges then fails with "Argument to filter condition must evaluate to -# a logical vector". - -#' @exportS3Method dplyr::filter -filter.RangedTupleList <- function(.data, ...) { - .rtlRequirePlyranges("filter") - out <- dplyr::filter(flattenTupleRanges(.data), ...) - nestTupleRanges(.rtlAsRanges(out), .data) -} - -#' @exportS3Method dplyr::mutate -mutate.RangedTupleList <- function(.data, ...) { - .rtlRequirePlyranges("mutate") - out <- dplyr::mutate(flattenTupleRanges(.data), ...) - nestTupleRanges(.rtlAsRanges(out), .data) -} - -#' @exportS3Method dplyr::select -select.RangedTupleList <- function(.data, ...) { - .rtlRequirePlyranges("select") - # The identity columns are kept whatever the selection, or the result - # could not be nested back into its elements. - out <- dplyr::select(flattenTupleRanges(.data), ...) - nestTupleRanges( - .rtlRestoreKeys(.rtlAsRanges(out), flattenTupleRanges(.data), .data), - .data - ) -} - -# Takes n ranges FROM EACH element, matching arrange()'s within-element rule. -# -# Unlike the other verbs this does not flatten: slicing the flat set would take -# n ranges in TOTAL, which for a many-element collection empties all but the -# first. Applying the slice per element is both the right semantics and simpler -# than reconstructing the grouping after a flatten. -#' @exportS3Method dplyr::slice -slice.RangedTupleList <- function(.data, ...) { - .rtlRequirePlyranges("slice") - elements <- map(as.list(.data), .rtlSliceOne, ...) - .rtlRebuild(.data, elements, seq_len(length(.data))) -} - -# @noRd -.rtlSliceOne <- function(g, ...) { - dplyr::slice(g, ...) -} - -# Orders WITHIN each element. A global sort of the flattened set followed by -# partitioning on the identity tuple leaves each element internally sorted -- -# the two are the same thing for the partitioned result -- so no grouping is -# needed here. Ordering ACROSS elements would permute the elements themselves, -# which the identity tuple cannot express. -#' @exportS3Method dplyr::arrange -arrange.RangedTupleList <- function(.data, ...) { - .rtlRequirePlyranges("arrange") - out <- dplyr::arrange(flattenTupleRanges(.data), ...) - nestTupleRanges(.rtlAsRanges(out), .data) -} - -#' @exportS3Method dplyr::summarise -summarise.RangedTupleList <- function(.data, ...) { - .rtlRequirePlyranges("summarise") - .rtlAsRanges(dplyr::summarise(flattenTupleRanges(.data), ...)) -} - -# No count() method: plyranges defines none for GRanges either, so there is -# nothing to delegate to. Use summarise(group_by(x, ...), n = n()). - -#' @exportS3Method dplyr::group_by -group_by.RangedTupleList <- function(.data, ...) { - .rtlRequirePlyranges("group_by") - .rtlAsRanges(dplyr::group_by(flattenTupleRanges(.data), ...)) -} - -# Put back any identity column a selection dropped. Without them the ranges -# carry no tuple and every one would fall into the first element. -# @noRd -.rtlRestoreKeys <- function(out, flat, template) { - dotted <- map_chr(.rtlTupleKeyCols(template), .rtlDotName) - missing <- setdiff( - intersect(dotted, colnames(mcols(flat))), - colnames(mcols(out)) - ) - for (nm in missing) { - mcols(out)[[nm]] <- mcols(flat)[[nm]] - } - out -} diff --git a/R/fineMappingPipeline.R b/R/fineMappingPipeline.R index aa498b3f..15fc4450 100644 --- a/R/fineMappingPipeline.R +++ b/R/fineMappingPipeline.R @@ -1027,18 +1027,10 @@ setGeneric("fineMappingPipeline", function(data, ...) { ".rbindFineMappingResult expects two FineMappingResultBase inputs." ) } - if (!identical(class(a)[[1L]], class(b)[[1L]])) { - clsA <- class(a)[[1L]] - clsB <- class(b)[[1L]] - msg <- glue( - ".rbindFineMappingResult: inputs must be the same concrete ", - "class (got '{clsA}' and '{clsB}')." - ) - abort(msg) - } - # Carry forward every column (blockId / joint* / ...) via the generic - # combine; the concrete class (QTL vs GWAS) is preserved automatically. - .rbindCollections(list(a, b), ldSketch = ldSketch) + # Carry forward every column (blockId / joint* / ...) and reconcile the + # collection-level slots via the shared combine; the concrete class + # (QTL vs GWAS) is preserved and checked there. + .combineTupleCollections(list(a, b), ldSketch, ".rbindFineMappingResult") } #' Combine FineMappingResult collections @@ -1067,7 +1059,7 @@ combineFineMappingResults <- function(..., ldSketch = NULL) { "FineMappingResultBase", "combineFineMappingResults" ) - reduce(parts, .rbindFineMappingResult, ldSketch = ldSketch) + .combineTupleCollections(parts, ldSketch, "combineFineMappingResults") } diff --git a/R/fineMappingWrappers.R b/R/fineMappingWrappers.R index a0ee4506..bed515f6 100644 --- a/R/fineMappingWrappers.R +++ b/R/fineMappingWrappers.R @@ -3310,6 +3310,245 @@ mergeSusieCs <- function(fineMappingResult, coverage = 0.95) { } +# ============================================================================= +# Post-fine-mapping credible-set extraction +# ----------------------------------------------------------------------------- +# Pull the trimmed SuSiE fit out of a pipeline result and reduce it to the +# per-credible-set / top-PIP diagnostic rows that summaryStatsQc() and the +# fine-mapping report consume. Relocated here from sumstatsQc.R: these read a +# fine-mapping result, they do not perform QC. +# ============================================================================= + +#' Extract the trimmed SuSiE fit from a finemapping pipeline result +#' +#' Returns the trimmed model fit underlying \code{con_data$finemappingEntry} (a +#' \code{FineMappingRow} S4 object), or NULL if no fine-mapping entry is +#' attached. +#' +#' @param conData List. The method-layer entry from a finemapping pipeline +#' result, expected to carry \code{$finemappingEntry} as a +#' \code{FineMappingRow} object. +#' @return The trimmed fit (a list with \code{pip}, \code{sets}, etc.) or NULL. +#' @examples +#' data(qtlSumStatsExample) +#' getSusieResult(qtlSumStatsExample) +#' @export +getSusieResult <- function(conData) { + if (length(conData) == 0) { + return(NULL) + } + fm <- conData$finemappingEntry + if (is.null(fm) || !is(fm, "FineMappingResultBase")) { + return(NULL) + } + trimmed <- .fmrPartsSusieFit(fm) + if (length(trimmed) == 0) { + return(NULL) + } + trimmed +} + +#' Process Credible Sets (CS) from Finemapping Results +#' +#' This function extracts and processes information for each Credible Set (CS) +#' from finemapping results, typically obtained from a finemapping RDS file. +#' +#' @param fmRow A \code{\link{fineMappingRow}} carrying the SuSiE +#' fit and variant ids (e.g. from \code{\link{getFineMappingResult}}). +#' @param csNames Character vector. Names of the Credible Sets, usually in the +#' format "L_". +#' @param topLociTable Data frame. The top-loci table (e.g. from +#' \code{\link{getTopLoci}}) carrying \code{variant_id}, \code{pip}, and +#' \code{z} columns. +#' @param ldSource The LD source from which the between-credible-set correlation +#' is derived on demand: a \code{QtlDataset} (individual-level) or a +#' \code{QtlSumStats} / \code{GwasSumStats} (summary statistics). See +#' \code{\link{computeCsCorrelation}}. +#' +#' @return A data frame with one row per CS, containing the following columns: +#' \item{cs_name}{Name of the Credible Set} +#' \item{variants_per_cs}{Number of variants in the CS} +#' \item{top_variant}{ID of the variant with the highest PIP in the CS} +#' \item{top_variant_index}{Global index of the top variant} +#' \item{top_pip}{Highest Posterior Inclusion Probability (PIP) in the CS} +#' \item{top_z}{Z-score of the top variant} +#' \item{p_value}{P-value calculated from the top Z-score} +#' \item{cs_corr_1, cs_corr_2, ...}{Each CS's pairwise correlation with every +#' CS (its row of the between-CS matrix, self-correlation on the diagonal), +#' computed on demand from \code{ldSource}. Absent when there are fewer than +#' two credible sets.} +#' \item{cs_corr_max}{Maximum absolute between-CS correlation (excluding the +#' self == 1); \code{NA} for a single CS.} +#' \item{cs_corr_min}{Minimum absolute between-CS correlation; \code{NA} for a +#' single CS.} +#' +#' @details This function is designed to be used only when there is at least one +#' Credible Set in the finemapping results usually for a given study and +#' block. It processes each CS, extracting key information such as the top +#' variant, its statistics, and correlation information between multiple CS if +#' available. +#' +#' @importFrom purrr map map_dbl map_int +#' @importFrom dplyr bind_rows +#' +#' @examples +#' data(qtlSumStatsExample) +#' vids <- c("chr1:100:A:G", "chr1:200:C:T", "chr1:300:G:A") +#' fit <- list(pip = c(0.1, 0.7, 0.2), sets = list(cs = list(L_1 = c(1, 2)))) +#' tl <- data.frame(variant_id = vids, pip = c(0.1, 0.7, 0.2), +#' z = c(1.0, 3.5, -0.5)) +#' fe <- fineMappingRow(variantIds = vids, susieFit = fit, topLoci = tl) +#' # A single credible set has no between-CS correlation (cs_corr_* are NA), so +#' # the ldSource is not consulted here. +#' extractCsInfo(fe, csNames = "L_1", topLociTable = tl, +#' ldSource = qtlSumStatsExample) +#' +#' @export +extractCsInfo <- function(fmRow, csNames, topLociTable, ldSource) { + fm <- fmRow + trimmed <- .fmrPartsSusieFit(fm) + variantNames <- .fmrPartsVariantIds(fm) + csCorr <- .rowCsCorrelation(fm, ldSource) + rows <- map( + seq_along(csNames), + .extractCsInfoRow, + csNames = csNames, + trimmed = trimmed, + variantNames = variantNames, + topLociTable = topLociTable + ) + .csAppendCorrelationCols(bind_rows(rows), csCorr) +} + +#' Extract Information for Top Variant from Finemapping Results +#' +#' This function extracts information about the variant with the highest +#' Posterior Inclusion Probability (PIP) from finemapping results, typically +#' used when no Credible Sets (CS) are identified in the analysis. +#' +#' @param fmRow A \code{\link{fineMappingRow}} carrying the SuSiE +#' fit and variant ids (e.g. from \code{\link{getFineMappingResult}}). +#' @param sumstats A list or data frame carrying a \code{z} element aligned to +#' the fit's variants (\code{sumstats$z}). +#' +#' @return A data frame with one row containing the following columns: +#' \item{cs_name}{NA (as no CS is identified)} +#' \item{variants_per_cs}{NA (as no CS is identified)} +#' \item{top_variant}{ID of the variant with the highest PIP} +#' \item{top_variant_index}{Index of the top variant in the original data} +#' \item{top_pip}{Highest Posterior Inclusion Probability (PIP)} +#' \item{top_z}{Z-score of the top variant} +#' \item{p_value}{P-value calculated from the top Z-score} +#' \item{cs_corr_max}{NA (no between-CS correlation without a CS)} +#' \item{cs_corr_min}{NA (no between-CS correlation without a CS)} +#' +#' @details This function is designed to be used when no Credible Sets are +#' identified in the finemapping results, but information about the most +#' significant variant is still desired. It identifies the variant with the +#' highest PIP and extracts relevant statistical information. +#' +#' @note This function is particularly useful for capturing information about +#' potentially important variants that might be included in Credible Sets +#' under different analysis parameters or lower coverage. It maintains a +#' structure similar to the output of `extract_cs_info()` for consistency in +#' downstream analyses. +#' +#' @seealso \code{\link{extractCsInfo}} for processing when Credible Sets are +#' present. +#' +#' @examples +#' vids <- c("chr1:100:A:G", "chr1:200:C:T", "chr1:300:G:A") +#' fit <- list(pip = c(0.1, 0.7, 0.2)) +#' tl <- data.frame(variant_id = vids, pip = c(0.1, 0.7, 0.2)) +#' fe <- fineMappingRow(variantIds = vids, susieFit = fit, topLoci = tl) +#' extractTopPipInfo(fe, sumstats = list(z = c(1.0, 3.5, -0.5))) +#' +#' @export +extractTopPipInfo <- function(fmRow, sumstats) { + fm <- fmRow + trimmed <- .fmrPartsSusieFit(fm) + variantNames <- .fmrPartsVariantIds(fm) + # Find the variant with the highest PIP + topPipIndex <- which.max(trimmed$pip) + topPip <- trimmed$pip[topPipIndex] + topVariant <- variantNames[topPipIndex] + topZ <- sumstats$z[topPipIndex] + pValue <- .zToPvalue(topZ) + + list( + cs_name = NA, + variants_per_cs = NA, + top_variant = topVariant, + top_variant_index = topPipIndex, + top_pip = topPip, + top_z = topZ, + p_value = pValue, + cs_corr_max = NA_real_, + cs_corr_min = NA_real_ + ) +} + +# Reduce one credible set's correlation vector (a row of the between-CS matrix, +# the self-correlation == 1 on the diagonal) to its |corr| max/min, excluding +# every self / perfect correlation (== 1). An empty result yields NA. +# @noRd +.extractCorrelations <- function(x) { + filtered <- abs(x[x != 1]) + if (length(filtered) == 0L) { + return(list(max_corr = NA_real_, min_corr = NA_real_)) + } + list( + max_corr = max(filtered, na.rm = TRUE), + min_corr = min(filtered, na.rm = TRUE) + ) +} + +# Append the between-CS correlation columns to the per-CS summary `base` from +# the m x m matrix `csCorr` (whose rows are aligned to `base`): the expanded +# cs_corr_1..m (each CS's row, self-correlation on the diagonal) plus +# cs_corr_max / cs_corr_min (|corr| excluding the self == 1). A NULL matrix +# (fewer than two credible sets) yields NA max/min and no expanded columns. +# @noRd +.csAppendCorrelationCols <- function(base, csCorr) { + if (is.null(csCorr)) { + return(mutate(base, cs_corr_max = NA_real_, cs_corr_min = NA_real_)) + } + perRow <- apply(csCorr, 1L, .extractCorrelations, simplify = FALSE) + expanded <- as_tibble(csCorr, .name_repair = "minimal") + names(expanded) <- str_c("cs_corr_", seq_len(ncol(csCorr))) + # unname(): apply() names its result by the matrix rownames, which map_dbl + # then carries into the column (tibbles preserve element names). + base |> + bind_cols(expanded) |> + mutate( + cs_corr_max = unname(map_dbl(perRow, "max_corr")), + cs_corr_min = unname(map_dbl(perRow, "min_corr")) + ) +} + +# One credible set's scalar summary row (top variant / PIP / z / p). The +# between-CS correlation columns are appended once by .csAppendCorrelationCols. +# @noRd +.extractCsInfoRow <- function(i, csNames, trimmed, variantNames, topLociTable) { + csName <- csNames[i] + indices <- trimmed$sets$cs[[csName]] + csVariants <- variantNames[indices] + csData <- filter(topLociTable, is_in(.data$variant_id, csVariants)) + topRow <- which.max(csData$pip) + topVariant <- csData$variant_id[topRow] + topZ <- csData$z[topRow] + tibble( + cs_name = csName, + variants_per_cs = length(csVariants), + top_variant = topVariant, + top_variant_index = which(variantNames == topVariant), + top_pip = csData$pip[topRow], + top_z = topZ, + p_value = .zToPvalue(topZ) + ) +} + + # ============================================================================= # SuSiE-family fitters (single-fit wrappers + per-block dispatch) # ----------------------------------------------------------------------------- diff --git a/R/gwasSumStats.R b/R/gwasSumStats.R index c5f25b53..876f039b 100644 --- a/R/gwasSumStats.R +++ b/R/gwasSumStats.R @@ -11,6 +11,18 @@ #' @include AllClasses.R tupleSelectors.R NULL +#' @title GWAS Summary-Statistic Collection +#' @description S4 collection of GWAS summary statistics keyed by the identity +#' tuple \code{(study)}. Each element is that study's per-variant +#' \code{GRanges} covering a single LD block, so build one collection per +#' block when sweeping the genome. +#' @details Required column: \code{study}, unique across rows. The class-level +#' slots inherited from \code{\linkS4class{SumStatsBase}} -- +#' \code{ldSketch}, \code{genome} and \code{qcInfo} -- apply uniformly to +#' every row. +#' @seealso \code{\link{GwasSumStats}} for the constructor and +#' \code{\linkS4class{QtlSumStats}} for the QTL counterpart. +#' @export setClass( "GwasSumStats", contains = "SumStatsBase", @@ -25,12 +37,7 @@ setClass( str_c("missing columns: ", str_flatten(missingCols, ", ")) ) } - if (length(object@genome) != 1L || str_length(object@genome) == 0L) { - errors <- c( - errors, - "'genome' slot must be a single non-empty character string" - ) - } + errors <- c(errors, .sumStatsCheckGenome(object)) if (!is.list(object@qcInfo)) { errors <- c(errors, "'qcInfo' slot must be a list") } @@ -58,7 +65,7 @@ setClass( setMethod("show", "GwasSumStats", function(object) { cat(glue( "GwasSumStats: {nrow(object)} studies, ", - "genome build {object@genome}\n", + "genome build {getGenome(object)}\n", .trim = FALSE )) ld <- object@ldSketch @@ -215,11 +222,14 @@ GwasSumStats <- function( # coerce those slots. Used by qtlSumStats.R too. # @noRd .sumStatsNewValidated <- function(Class, grl, ldSketch, genome, qcInfo) { + # The build goes into seqinfo, which is where a GRangesList keeps it and + # where every Bioconductor consumer reads it from. Assigned before new() + # so validity sees the finished object. + GenomeInfoDb::genome(grl) <- as.character(genome) obj <- methods::new( Class, grl, ldSketch = .asLdSketch(ldSketch), - genome = as.character(genome), qcInfo = as.list(qcInfo) ) methods::validObject(obj) diff --git a/R/qtlAssociationPostprocess.R b/R/qtlAssociationPostprocess.R index 12059f02..db9778f6 100644 --- a/R/qtlAssociationPostprocess.R +++ b/R/qtlAssociationPostprocess.R @@ -283,12 +283,15 @@ setMethod( md[[nm]] <- newCols[[nm]] } grl <- GenomicRanges::GRangesList(as.list(x)) + # Rebuilding from as.list() starts from the elements' own seqinfo, so the + # build is written back explicitly -- it is collection-level state, and + # there is no genome slot to carry it any more. + GenomeInfoDb::genome(grl) <- getGenome(x) mcols(grl) <- md methods::new( "QtlSumStats", grl, ldSketch = getLdSketch(x), - genome = getGenome(x), qcInfo = qcInfo ) } diff --git a/R/qtlSumStats.R b/R/qtlSumStats.R index 2dd4026e..488c4974 100644 --- a/R/qtlSumStats.R +++ b/R/qtlSumStats.R @@ -12,6 +12,19 @@ #' @include AllClasses.R tupleSelectors.R NULL +#' @title QTL Summary-Statistic Collection +#' @description S4 collection of QTL summary statistics keyed by the identity +#' tuple \code{(study, context, trait)}. Each element is that tuple's +#' per-variant \code{GRanges} -- \code{variant_id} plus the per-variant +#' Z / N / MAF mcols. +#' @details Required columns: \code{study}, \code{context} and \code{trait}; +#' the 3-tuple is unique. The class-level slots inherited from +#' \code{\linkS4class{SumStatsBase}} -- \code{ldSketch}, \code{genome} and +#' \code{qcInfo} -- apply uniformly to every row, and \code{qcInfo} records +#' which \code{\link{summaryStatsQc}} passes have run. +#' @seealso \code{\link{QtlSumStats}} for the constructor and +#' \code{\linkS4class{GwasSumStats}} for the GWAS counterpart. +#' @export setClass( "QtlSumStats", contains = "SumStatsBase", @@ -58,13 +71,10 @@ setClass( NULL } -# genome slot: a single non-empty string. +# The genome build, read from seqinfo (there is no genome slot). # @noRd .qssCheckGenome <- function(object) { - if (length(object@genome) != 1L || str_length(object@genome) == 0L) { - return("'genome' slot must be a single non-empty character string") - } - NULL + .sumStatsCheckGenome(object) } # qcInfo slot must be a list. diff --git a/R/sumstatsQc.R b/R/sumstatsQc.R index 5bd26c6c..1e907463 100644 --- a/R/sumstatsQc.R +++ b/R/sumstatsQc.R @@ -1390,216 +1390,9 @@ slalom <- function( # ============================================================================= -# Univariate RSS diagnostics (post-finemap) +# Post-finemap credible-set update strategy # ============================================================================= -#' Extract the trimmed SuSiE fit from a finemapping pipeline result -#' -#' Returns the trimmed model fit underlying \code{con_data$finemappingEntry} (a -#' \code{FineMappingRow} S4 object), or NULL if no fine-mapping entry is -#' attached. -#' -#' @param conData List. The method-layer entry from a finemapping pipeline -#' result, expected to carry \code{$finemappingEntry} as a -#' \code{FineMappingRow} object. -#' @return The trimmed fit (a list with \code{pip}, \code{sets}, etc.) or NULL. -#' @examples -#' data(qtlSumStatsExample) -#' getSusieResult(qtlSumStatsExample) -#' @export -getSusieResult <- function(conData) { - if (length(conData) == 0) { - return(NULL) - } - fm <- conData$finemappingEntry - if (is.null(fm) || !is(fm, "FineMappingResultBase")) { - return(NULL) - } - trimmed <- .fmrPartsSusieFit(fm) - if (length(trimmed) == 0) { - return(NULL) - } - trimmed -} - -#' Process Credible Sets (CS) from Finemapping Results -#' -#' This function extracts and processes information for each Credible Set (CS) -#' from finemapping results, typically obtained from a finemapping RDS file. -#' -#' @param fmRow A \code{\link{fineMappingRow}} carrying the SuSiE -#' fit and variant ids (e.g. from \code{\link{getFineMappingResult}}). -#' @param csNames Character vector. Names of the Credible Sets, usually in the -#' format "L_". -#' @param topLociTable Data frame. The top-loci table (e.g. from -#' \code{\link{getTopLoci}}) carrying \code{variant_id}, \code{pip}, and -#' \code{z} columns. -#' @param ldSource The LD source from which the between-credible-set correlation -#' is derived on demand: a \code{QtlDataset} (individual-level) or a -#' \code{QtlSumStats} / \code{GwasSumStats} (summary statistics). See -#' \code{\link{computeCsCorrelation}}. -#' -#' @return A data frame with one row per CS, containing the following columns: -#' \item{cs_name}{Name of the Credible Set} -#' \item{variants_per_cs}{Number of variants in the CS} -#' \item{top_variant}{ID of the variant with the highest PIP in the CS} -#' \item{top_variant_index}{Global index of the top variant} -#' \item{top_pip}{Highest Posterior Inclusion Probability (PIP) in the CS} -#' \item{top_z}{Z-score of the top variant} -#' \item{p_value}{P-value calculated from the top Z-score} -#' \item{cs_corr_1, cs_corr_2, ...}{Each CS's pairwise correlation with every -#' CS (its row of the between-CS matrix, self-correlation on the diagonal), -#' computed on demand from \code{ldSource}. Absent when there are fewer than -#' two credible sets.} -#' \item{cs_corr_max}{Maximum absolute between-CS correlation (excluding the -#' self == 1); \code{NA} for a single CS.} -#' \item{cs_corr_min}{Minimum absolute between-CS correlation; \code{NA} for a -#' single CS.} -#' -#' @details This function is designed to be used only when there is at least one -#' Credible Set in the finemapping results usually for a given study and -#' block. It processes each CS, extracting key information such as the top -#' variant, its statistics, and correlation information between multiple CS if -#' available. -#' -#' @importFrom purrr map map_dbl map_int -#' @importFrom dplyr bind_rows -#' -#' @examples -#' data(qtlSumStatsExample) -#' vids <- c("chr1:100:A:G", "chr1:200:C:T", "chr1:300:G:A") -#' fit <- list(pip = c(0.1, 0.7, 0.2), sets = list(cs = list(L_1 = c(1, 2)))) -#' tl <- data.frame(variant_id = vids, pip = c(0.1, 0.7, 0.2), -#' z = c(1.0, 3.5, -0.5)) -#' fe <- fineMappingRow(variantIds = vids, susieFit = fit, topLoci = tl) -#' # A single credible set has no between-CS correlation (cs_corr_* are NA), so -#' # the ldSource is not consulted here. -#' extractCsInfo(fe, csNames = "L_1", topLociTable = tl, -#' ldSource = qtlSumStatsExample) -#' -#' @export -extractCsInfo <- function(fmRow, csNames, topLociTable, ldSource) { - fm <- fmRow - trimmed <- .fmrPartsSusieFit(fm) - variantNames <- .fmrPartsVariantIds(fm) - csCorr <- .rowCsCorrelation(fm, ldSource) - rows <- map( - seq_along(csNames), - .extractCsInfoRow, - csNames = csNames, - trimmed = trimmed, - variantNames = variantNames, - topLociTable = topLociTable - ) - .csAppendCorrelationCols(bind_rows(rows), csCorr) -} - -#' Extract Information for Top Variant from Finemapping Results -#' -#' This function extracts information about the variant with the highest -#' Posterior Inclusion Probability (PIP) from finemapping results, typically -#' used when no Credible Sets (CS) are identified in the analysis. -#' -#' @param fmRow A \code{\link{fineMappingRow}} carrying the SuSiE -#' fit and variant ids (e.g. from \code{\link{getFineMappingResult}}). -#' @param sumstats A list or data frame carrying a \code{z} element aligned to -#' the fit's variants (\code{sumstats$z}). -#' -#' @return A data frame with one row containing the following columns: -#' \item{cs_name}{NA (as no CS is identified)} -#' \item{variants_per_cs}{NA (as no CS is identified)} -#' \item{top_variant}{ID of the variant with the highest PIP} -#' \item{top_variant_index}{Index of the top variant in the original data} -#' \item{top_pip}{Highest Posterior Inclusion Probability (PIP)} -#' \item{top_z}{Z-score of the top variant} -#' \item{p_value}{P-value calculated from the top Z-score} -#' \item{cs_corr_max}{NA (no between-CS correlation without a CS)} -#' \item{cs_corr_min}{NA (no between-CS correlation without a CS)} -#' -#' @details This function is designed to be used when no Credible Sets are -#' identified in the finemapping results, but information about the most -#' significant variant is still desired. It identifies the variant with the -#' highest PIP and extracts relevant statistical information. -#' -#' @note This function is particularly useful for capturing information about -#' potentially important variants that might be included in Credible Sets -#' under different analysis parameters or lower coverage. It maintains a -#' structure similar to the output of `extract_cs_info()` for consistency in -#' downstream analyses. -#' -#' @seealso \code{\link{extractCsInfo}} for processing when Credible Sets are -#' present. -#' -#' @examples -#' vids <- c("chr1:100:A:G", "chr1:200:C:T", "chr1:300:G:A") -#' fit <- list(pip = c(0.1, 0.7, 0.2)) -#' tl <- data.frame(variant_id = vids, pip = c(0.1, 0.7, 0.2)) -#' fe <- fineMappingRow(variantIds = vids, susieFit = fit, topLoci = tl) -#' extractTopPipInfo(fe, sumstats = list(z = c(1.0, 3.5, -0.5))) -#' -#' @export -extractTopPipInfo <- function(fmRow, sumstats) { - fm <- fmRow - trimmed <- .fmrPartsSusieFit(fm) - variantNames <- .fmrPartsVariantIds(fm) - # Find the variant with the highest PIP - topPipIndex <- which.max(trimmed$pip) - topPip <- trimmed$pip[topPipIndex] - topVariant <- variantNames[topPipIndex] - topZ <- sumstats$z[topPipIndex] - pValue <- .zToPvalue(topZ) - - list( - cs_name = NA, - variants_per_cs = NA, - top_variant = topVariant, - top_variant_index = topPipIndex, - top_pip = topPip, - top_z = topZ, - p_value = pValue, - cs_corr_max = NA_real_, - cs_corr_min = NA_real_ - ) -} - -# Reduce one credible set's correlation vector (a row of the between-CS matrix, -# the self-correlation == 1 on the diagonal) to its |corr| max/min, excluding -# every self / perfect correlation (== 1). An empty result yields NA. -# @noRd -.extractCorrelations <- function(x) { - filtered <- abs(x[x != 1]) - if (length(filtered) == 0L) { - return(list(max_corr = NA_real_, min_corr = NA_real_)) - } - list( - max_corr = max(filtered, na.rm = TRUE), - min_corr = min(filtered, na.rm = TRUE) - ) -} - -# Append the between-CS correlation columns to the per-CS summary `base` from -# the m x m matrix `csCorr` (whose rows are aligned to `base`): the expanded -# cs_corr_1..m (each CS's row, self-correlation on the diagonal) plus -# cs_corr_max / cs_corr_min (|corr| excluding the self == 1). A NULL matrix -# (fewer than two credible sets) yields NA max/min and no expanded columns. -# @noRd -.csAppendCorrelationCols <- function(base, csCorr) { - if (is.null(csCorr)) { - return(mutate(base, cs_corr_max = NA_real_, cs_corr_min = NA_real_)) - } - perRow <- apply(csCorr, 1L, .extractCorrelations, simplify = FALSE) - expanded <- as_tibble(csCorr, .name_repair = "minimal") - names(expanded) <- str_c("cs_corr_", seq_len(ncol(csCorr))) - # unname(): apply() names its result by the matrix rownames, which map_dbl - # then carries into the column (tibbles preserve element names). - base |> - bind_cols(expanded) |> - mutate( - cs_corr_max = unname(map_dbl(perRow, "max_corr")), - cs_corr_min = unname(map_dbl(perRow, "min_corr")) - ) -} - #' Process Credible Set Information and Determine Updating Strategy #' #' This function categorizes Credible Sets (CS) within a study block into @@ -5197,28 +4990,6 @@ summaryStatsQc <- function( # ---- other map/apply helpers (lambda-free callbacks) -------------------- -# One credible set's scalar summary row (top variant / PIP / z / p). The -# between-CS correlation columns are appended once by .csAppendCorrelationCols. -# @noRd -.extractCsInfoRow <- function(i, csNames, trimmed, variantNames, topLociTable) { - csName <- csNames[i] - indices <- trimmed$sets$cs[[csName]] - csVariants <- variantNames[indices] - csData <- filter(topLociTable, is_in(.data$variant_id, csVariants)) - topRow <- which.max(csData$pip) - topVariant <- csData$variant_id[topRow] - topZ <- csData$z[topRow] - tibble( - cs_name = csName, - variants_per_cs = length(csVariants), - top_variant = topVariant, - top_variant_index = which(variantNames == topVariant), - top_pip = csData$pip[topRow], - top_z = topZ, - p_value = .zToPvalue(topZ) - ) -} - # TRUE when credible set row `i` is a tagged (redundant) set. # @noRd .autoDecisionTagged <- function(i, df, highCorrCols) { diff --git a/R/tupleSelectors.R b/R/tupleSelectors.R index 8fa0cbb1..436699cc 100644 --- a/R/tupleSelectors.R +++ b/R/tupleSelectors.R @@ -1,9 +1,11 @@ # ============================================================================= # Tuple row matchers # ----------------------------------------------------------------------------- -# Internal helpers shared by the FineMappingResult / TwasWeights / -# SumStats DFrame-subclass collections to resolve a tuple-keyed selection -# to a single row index. Pure R helpers -- no S4 dispatch, no exports. +# Internal helpers shared by the FineMappingResult / TwasWeights / SumStats +# collections to resolve a tuple-keyed selection to a single row index. Every +# collection is a RangedTupleList, so these read identity columns off mcols +# and payloads off the elements -- there is no second shape to branch on. +# Pure R helpers -- no S4 dispatch, no exports. # ============================================================================= # Internal: return integer row indices of `x` where every (column, value) @@ -23,16 +25,12 @@ which(ok) } -# Read an identity column by name, whichever collection shape `x` has. -# `x[[k]]` is a column on the DFrame-backed collections but an ELEMENT on a -# RangedTupleList, where the identity columns live in mcols. The two shapes -# coexist until every collection has migrated. +# Read an identity column by name. Every collection keeps its identity +# columns in mcols; `x[[k]]` addresses an ELEMENT, not a column, so mcols is +# the only correct read. # @noRd .tupleColumn <- function(x, k) { - if (methods::is(x, "RangedTupleList")) { - return(mcols(x)[[k]]) - } - x[[k]] + mcols(x)[[k]] } # One collection ELEMENT, whichever shape `x` has. On a RangedTupleList the @@ -47,10 +45,7 @@ if (methods::is(x, "FineMappingResultBase")) { return(.fmrRowParts(x, i)) } - if (methods::is(x, "RangedTupleList")) { - return(x[[i]]) - } - x$entry[[i]] + x[[i]] } # ALL elements as a plain list, whichever shape `x` has. @@ -64,20 +59,14 @@ ) { return(map(seq_len(nrow(x)), .collectionEntry, x = x)) } - if (methods::is(x, "RangedTupleList")) { - return(as.list(x)) - } - as.list(x$entry) + as.list(x) } -# The identity column NAMES, whichever shape `x` has. `names(x)` is element -# names on a RangedTupleList, not columns. +# The identity column NAMES. `names(x)` would be the ELEMENT names, not the +# columns, so this reads mcols. # @noRd .tupleColumnNames <- function(x) { - if (methods::is(x, "RangedTupleList")) { - return(colnames(mcols(x))) - } - names(x) + colnames(mcols(x)) } # Internal: resolve a tuple-keyed selection (study, context, trait, @@ -286,13 +275,9 @@ if (nrow(x) == 0L) { return(GenomicRanges::GRanges()) } - if (methods::is(x, "RangedTupleList")) { - return(unlist(range(x), use.names = FALSE)) - } - # No fallback to a stored column: nothing carries one any more. A DFrame - # collection (CtwasResult) has no region column either, so an empty - # GRanges is the honest answer rather than a fabricated span. - GenomicRanges::GRanges() + # No stored region column: nothing carries one any more, so the span is + # always derived from the elements' own ranges. + unlist(range(x), use.names = FALSE) } # Internal: append the optional `blockId` provenance column (no-op when NULL). @@ -406,30 +391,26 @@ # region) is padded with `.naLikeColumn` -- so adding a column to a class flows # through combines automatically rather than being hand-listed at each site. # Preserves the concrete class and sets the `ldSketch` slot explicitly (rbind -# does not reliably carry it). Returns NULL when given no inputs. +# does not reliably carry it). Returns NULL when given no inputs. Every +# collection is a RangedTupleList now, so this only drops the NULL parts and +# hands off; the old DFrame path went with the rebase onto GRangesList. # # `slots` carries the collection-level state a concrete subclass adds beyond # `ldSketch` (SumStatsBase's `genome` / `qcInfo`). It is passed in rather than # read off `parts[[1L]]` because those slots describe the WHOLE collection, so # merging them is the caller's decision -- first-wins would silently drop the # other parts' QC audit. -.rbindCollections <- function(parts, ldSketch = NULL, slots = list()) { +.rbindCollections <- function( + parts, + ldSketch = NULL, + slots = list(), + genome = NULL +) { parts <- compact(parts) if (length(parts) == 0L) { return(NULL) } - if (methods::is(parts[[1L]], "RangedTupleList")) { - return(.rbindRangedCollections(parts, ldSketch, slots)) - } - cls <- class(parts[[1L]])[[1L]] - allCols <- reduce(map(parts, names), union) - combined <- map(allCols, .rbindColumn, parts = parts) - names(combined) <- allCols - dfArgs <- c(combined, list(check.names = FALSE)) - df <- exec(S4Vectors::DataFrame, !!!dfArgs) - out <- new(cls, df, ldSketch = ldSketch) - validObject(out) - out + .rbindRangedCollections(parts, ldSketch, slots, genome) } # Internal: aggregate a per-entry accessor across every row of a @@ -588,11 +569,20 @@ # Concatenate column `cn` across all parts (NA-filling parts that lack it). # @noRd -# Row-binding a RANGED collection: the elements append and the mcols rbind. -# The DFrame path cannot be reused because there is no column holding the -# payload any more -- the payload is the container. +# Row-binding a collection: the elements append and the mcols rbind. +# +# S4Vectors::combineRows() does the union-with-NA-fill this needs, but it +# fails on a GRanges-valued mcols column (`traitPos`) with "GRanges objects +# don't support [[, as.list(), lapply(), or unlist()". .rbindColumn() binds +# each column with c() instead, which GRanges does support, so the union is +# done per column here rather than delegated. # @noRd -.rbindRangedCollections <- function(parts, ldSketch, slots = list()) { +.rbindRangedCollections <- function( + parts, + ldSketch, + slots = list(), + genome = NULL +) { cls <- class(parts[[1L]])[[1L]] allCols <- reduce(map(parts, .tupleColumnNames), union) combined <- set_names(map(allCols, .rbindColumn, parts = parts), allCols) @@ -602,6 +592,13 @@ ) elements <- list_flatten(map(parts, as.list)) grl <- GenomicRanges::GRangesList(elements) + # Rebuilding from as.list() merges each part's seqinfo, so a part whose + # elements never had a build set would leave the result carrying both the + # real build and NA. The agreed build is written back explicitly to keep + # seqinfo single-valued, which validity requires. + if (!is.null(genome)) { + GenomeInfoDb::genome(grl) <- genome + } mcols(grl) <- md out <- exec( methods::new, @@ -1006,21 +1003,84 @@ # first -- and cTWAS computes its full-panel LD from exactly that slot. # ============================================================================= -# Row-bind summary-statistics parts, merging the three collection-level slots. -# `ldSketch` (when non-NULL) overrides the unioned panel. +# The concrete class of `x`, for map_chr over a list of collections. # @noRd -.rbindSumStats <- function(parts, ldSketch, fn) { - sketch <- ldSketch %||% .ssCombineSketch(parts, fn) +.firstClass <- function(x) class(x)[[1L]] + +# Every part must be the same concrete class: the combined object is built as +# `parts[[1L]]`'s class, so a mixed set would silently coerce the rest. +# @noRd +.rtlRequireSameClass <- function(parts, fn) { + classes <- unique(map_chr(parts, .firstClass)) + if (length(classes) > 1L) { + msg <- glue( + "{fn}: inputs must be the same concrete class (got ", + "{str_flatten(classes, ', ')})." + ) + abort(msg) + } + invisible(NULL) +} + +# The collection-level slots BEYOND ldSketch that a combine has to merge, for +# whatever concrete class `parts` are. A slot describes the whole collection, +# so first-wins would silently drop the other parts' state; each gets its own +# rule and an unrecognised slot is an error rather than a quiet default, so +# adding one to a class cannot slip through a combine unnoticed. +# @noRd +.rtlExtraSlots <- function(parts, fn) { + own <- setdiff(.rtlOwnSlots(parts[[1L]]), "ldSketch") + out <- list() + if (is_in("qcInfo", own)) { + out$qcInfo <- .ssCombineQcInfo(parts, fn) + } + unknown <- setdiff(own, names(out)) + if (length(unknown) > 0L) { + cls <- class(parts[[1L]])[[1L]] + msg <- glue( + "{fn}: no merge rule for the collection-level slot(s) ", + "{str_flatten(unknown, ', ')} on class '{cls}'." + ) + abort(msg) + } + out +} + +# Merge `parts` into one collection, reconciling every collection-level slot. +# The single path behind `c()` / `append()` and the four combine*() wrappers. +# A non-NULL `ldSketch` overrides the unioned panel. +# @noRd +.combineTupleCollections <- function(parts, ldSketch, fn) { + .rtlRequireSameClass(parts, fn) + # Forced before the rebuild: GRangesList() merges seqinfo and rejects + # mismatched builds itself, with a message about sequence-level genomes + # that says nothing about which inputs disagreed. + genome <- .rtlCombineGenome(parts, fn) .rbindCollections( parts, - ldSketch = sketch, - slots = list( - genome = .ssCombineGenome(parts, fn), - qcInfo = .ssCombineQcInfo(parts, fn) - ) + ldSketch = ldSketch %||% .combineLdSketch(parts, fn), + slots = .rtlExtraSlots(parts, fn), + genome = genome ) } +# The single build every part must agree on, or NULL for a family that does +# not carry one. Read through getGenome() so it works off seqinfo. +# @noRd +.rtlCombineGenome <- function(parts, fn) { + if (!methods::is(parts[[1L]], "SumStatsBase")) { + return(NULL) + } + .ssCombineGenome(parts, fn) +} + +# Row-bind summary-statistics parts, merging the three collection-level slots. +# `ldSketch` (when non-NULL) overrides the unioned panel. +# @noRd +.rbindSumStats <- function(parts, ldSketch, fn) { + .combineTupleCollections(parts, ldSketch, fn) +} + # One genome build per collection (every entry shares the LD sketch), so the # parts must agree. # @noRd @@ -1104,7 +1164,7 @@ # collection has no panel); a mix is an error, because silently keeping the # panels that exist would leave elements harmonized against nothing. # @noRd -.ssCombineSketch <- function(parts, fn) { +.combineLdSketch <- function(parts, fn) { sketches <- map(parts, getLdSketch) present <- !map_lgl(sketches, is.null) if (!any(present)) { @@ -1122,9 +1182,35 @@ .ssCheckSameSource(handles, fn) handle <- handles[[1L]] handle@snpInfo <- .ssUnionSnpInfo(handles) + handle@chromPaths <- .ssUnionChromPaths(handles, fn) .asLdSketch(handle) } +# Union the per-chromosome shard paths across handles. Two sketches trimmed +# from one genoMeta panel to different chromosomes share path/format/samples +# and differ ONLY here, so unioning is what lets their snpInfo rows keep +# routing: a row reaches its file by its own CHR. One chromosome mapping to +# two different files means the parts are not the same panel after all. +# @noRd +.ssUnionChromPaths <- function(handles, fn) { + out <- character(0) + for (h in handles) { + cp <- h@chromPaths + for (ch in names(cp)) { + if (is_in(ch, names(out)) && !identical(out[[ch]], cp[[ch]])) { + msg <- glue( + "{fn}: chromosome '{ch}' maps to two different genotype ", + "files across the inputs, so their LD sketches cannot ", + "be unioned." + ) + abort(msg) + } + out[[ch]] <- cp[[ch]] + } + } + out +} + # Every panel must read from the same file(s): the union keeps the first # handle and widens its snpInfo, and each row's `fileIdx` addresses a position # in THAT file, so rows from another panel would read the wrong variants. @@ -1144,11 +1230,16 @@ # @noRd .ssSameGenotypeSource <- function(h, first) { + # chromPaths is deliberately NOT compared: two sketches trimmed from the + # same genoMeta panel to different chromosomes differ there and nowhere + # else, and .ssUnionChromPaths() merges them (erroring on a genuine + # conflict). Everything below must match for the union to mean anything -- + # snpInfo rows carry file POSITIONS, which only make sense against the + # same files and the same sample axis. identical(h@path, first@path) && identical(h@format, first@format) && identical(h@nSamples, first@nSamples) && - identical(h@sampleIds, first@sampleIds) && - identical(h@chromPaths, first@chromPaths) + identical(h@sampleIds, first@sampleIds) } # The parts' snpInfo rows, de-duplicated by variant id and returned in genomic diff --git a/R/twasWeightsPipeline.R b/R/twasWeightsPipeline.R index 3a185e53..d20819b9 100644 --- a/R/twasWeightsPipeline.R +++ b/R/twasWeightsPipeline.R @@ -10,8 +10,9 @@ if (!is(a, "TwasWeights") || !is(b, "TwasWeights")) { abort(".rbindTwasWeights expects two TwasWeights inputs.") } - # Carry forward every column (joint*, region, ...) via the generic combine. - .rbindCollections(list(a, b), ldSketch = ldSketch) + # Carry forward every column (joint*, region, ...) and reconcile the + # collection-level slots via the shared combine. + .combineTupleCollections(list(a, b), ldSketch, ".rbindTwasWeights") } # Normalize combine() varargs: accept either N objects or a single list of @@ -63,7 +64,7 @@ #' @export combineTwasWeights <- function(..., ldSketch = NULL) { parts <- .asCombineList(list(...), "TwasWeights", "combineTwasWeights") - reduce(parts, .rbindTwasWeights, ldSketch = ldSketch) + .combineTupleCollections(parts, ldSketch, "combineTwasWeights") } # --- Multi-region (jointRegions) helpers for the QtlDataset method ---------- diff --git a/_pkgdown.yml b/_pkgdown.yml index 1862f12c..907d04aa 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -92,15 +92,20 @@ reference: desc: > S4 class definitions. The `-class` topic documents the slots and validity constraints; the matching constructor topic (next - section) documents the user-facing factory function. + section) documents the user-facing factory function. Virtual + classes have no constructor, and neither does `H2Estimate` -- + `estimateH2()` builds it. contents: - AnnotationMatrix-class - ColocBoostResult-class - ColocResult-class - ColocResultBase-class + - CtwasResult-class + - CtwasResultEntry-class - FineMappingResultBase-class - FineMappingRow-class - GwasFineMappingResult-class + - GwasSumStats-class - H2Estimate-class - LdData-class - LdEigen-class @@ -110,6 +115,7 @@ reference: - MultiStudyQtlDataset-class - QtlDataset-class - QtlFineMappingResult-class + - QtlSumStats-class - RangedTupleList-class - SldscData-class - SumStatsBase-class @@ -118,15 +124,16 @@ reference: - title: "Class constructors" desc: > - User-facing constructors. For classes that have both a `-class` - topic and a constructor topic, prefer the constructor for + User-facing constructors -- one per class, plus the two + alternative `QtlSumStats` builders. Every class here also has a + `-class` topic in the section above: prefer the constructor for day-to-day use; the `-class` topic is the authoritative slot reference. contents: - AnnotationMatrix - - CtwasResult - ColocBoostResult - ColocResult + - CtwasResult - CtwasResultEntry - fineMappingRow - GwasFineMappingResult @@ -158,88 +165,160 @@ reference: - title: "Class methods" desc: > - Accessor and behaviour methods grouped by the class they're - defined on. Generics with implementations on multiple classes - appear under each. S4 pipeline dispatch methods - (`fineMappingPipeline`, `colocboostPipeline`, - `twasWeightsPipeline`) are listed under "Pipelines" further - down rather than repeated per class. + Accessor and behaviour methods, grouped by the class each method is + defined on. A generic with methods on several classes is listed under + every one of them, so each block is the complete callable surface for + that class (a subclass also inherits everything listed under its base). + S4 pipeline dispatch methods (`fineMappingPipeline`, + `colocboostPipeline`, `twasWeightsPipeline`) are listed under + "Pipelines" further down rather than repeated per class. contents: - show-methods - - subtitle: "AnnotationMatrix" - contents: - - getBaseline - - getCandidates - - getGenome - - subtitle: "RangedTupleList" contents: - # The virtual base every collection inherits: container and metadata - # methods, plus the region/variant verbs and the plyranges bridge. - - RangedTupleList-methods - - subsetRegion - - intersectVariants + # The virtual base of every tuple collection: container and metadata + # methods, the region verb, and the plyranges bridge. The families + # that inherit it -- fine-mapping results, colocalization results, + # summary statistics, TWAS weights -- follow. - flattenTupleRanges - nestTupleRanges + - RangedTupleList-methods + - subsetRegion + + - subtitle: "FineMappingResultBase" + contents: - computeCsCorrelation - fsusieAffectedRegions - fsusieCredibleBand - getCredibleSetSummary - getCs - getLbf + - getLdSketch - getMarginalEffects - - getPip + - getMethodNames + - getRegion + - getRetainedMass + - getStudy - getSusieFit - getTopLoci + - getTraitPosition - getVariantIds + - intersectVariants - resolveWeights - - - subtitle: "FineMappingResultBase" - contents: - - getLdSketch - - getMethodNames - - getRegion - - getStudy - writeSumstatsVcf - - subtitle: "CtwasResult / CtwasResultEntry" + - subtitle: "QtlFineMappingResult" contents: - - getCtwasParam - - getFinemap - - getSusieAlpha + - getContexts + - getCvResult + - getFineMappingResult + - getPip + - getTraits - subtitle: "GwasFineMappingResult" contents: - getContexts - getFineMappingResult + - getPip - getTraits - - subtitle: "GwasSumStats" + - subtitle: "FineMappingRow" + contents: + - getCvResult + - getSusieFit + - getVariantIds + + - subtitle: "ColocResultBase / ColocResult / ColocBoostResult" + contents: + - "as.data.frame,ColocBoostResult-method" + - "as.data.frame,ColocResult-method" + - colocViews + - getColocBoostOutcomes + - getComputingTime + - getLdSketch + - getRegionVcp + + - subtitle: "SumStatsBase" contents: - - "getSumStats,GwasSumStats-method" - - combineGwasSumStats - getBeta + - getGenome + - getLdSketch - getMaf - getN - getP + - getQcDiagnostics + - getQcInfo - getSe - - getSumStats - - getSumstatDf - - getVarY + - getStudy + - getVariantIds - getZ - nSnps - subsetChr - - subtitle: "ColocResult / ColocBoostResult" + - subtitle: "GwasSumStats" contents: - - colocViews - - getColocBoostOutcomes - - getComputingTime - - getRegionVcp - - getRetainedMass - - "as.data.frame,ColocResult-method" - - "as.data.frame,ColocBoostResult-method" + - combineGwasSumStats + - getSumstatDf + - "getSumStats,GwasSumStats-method" + - getVarY + - writeSumstatsVcf + + - subtitle: "QtlSumStats" + contents: + - combineQtlSumStats + - getContexts + - getSignificantQtls + - getSumstatDf + - getSumStats + - getTraitPosition + - getTraits + - getVarY + - qtlAssociationPostprocess + + - subtitle: "TwasWeights" + contents: + - getContexts + - getCvResult + - getDataType + - getFits + - getLdSketch + - getMethodNames + - getRegion + - getStandardized + - getStudy + - getTraitPosition + - getTraits + - getTwasWeights + - getVariantIds + - getWeights + - resolveWeights + + - subtitle: "TwasWeightsRow" + contents: + - getCvResult + - getDataType + - getFits + - getStandardized + - getVariantIds + - getWeights + + - subtitle: "AnnotationMatrix" + contents: + # Classes outside the RangedTupleList hierarchy follow, in + # alphabetical order. + - getBaseline + - getCandidates + - getGenome + + - subtitle: "CtwasResult / CtwasResultEntry" + contents: + - getContexts + - getCtwasParam + - getFinemap + - getMethodNames + - getStudy + - getSusieAlpha - subtitle: "H2Estimate" contents: @@ -249,6 +328,7 @@ reference: - getIntercept - getInterceptSe - getLocal + - getMethodNames - getNSnps - getScoreStats - getTauBlocks @@ -263,6 +343,7 @@ reference: - getNRef - getRefPanel - getSnpIdx + - getVariantIds - getVariantInfo - hasGenotypes @@ -273,80 +354,51 @@ reference: - subtitle: "LdScore" contents: - getLdMatrixList - - getLdScoreWeights - getLdScores + - getLdScoreWeights - subtitle: "LdStatistic" contents: + - getGenome - getInSample - getLdBlocks + - getNRef + + - subtitle: "MashPrior" + contents: + - getCvFits + - getFullFit - subtitle: "MultiStudyQtlDataset" contents: - getQtlDatasets + - getStudy + - getSumStats - subtitle: "QtlDataset" contents: - getAf + - getContexts - getGenotypeCovariates - getGenotypes + - getMaf - getPhenotypeCovariates - getPhenotypes - getResidualizedGenotypes - getResidualizedPhenotypes - getScaleResiduals + - getStudy - getTraitPosition - qtlDatasetFilters - - subtitle: "QtlFineMappingResult" - contents: - - getFineMappingResult - - - subtitle: "QtlSumStats" - contents: - - combineQtlSumStats - - getContexts - - getSignificantQtls - - getTraits - - qtlAssociationPostprocess - - - subtitle: "SumStatsBase" - contents: - - getQcInfo - - getQcDiagnostics - - subtitle: "SldscData" contents: + - getAnnotCols - getAnnotData - getFrqData - - getTraitRuns - getTraitNames - - getAnnotCols - getTraitRun - - - subtitle: "MashPrior" - contents: - - getCvFits - - getFullFit - - - subtitle: "TwasWeights" - contents: - - getCvResult - - getDataType - - getFits - - getStandardized - - getTwasWeights - - getWeights - - - subtitle: "TwasWeightsEntry" - contents: - # Same generic surface as TwasWeights — dispatches on per-row entries. - - getCvResult - - getDataType - - getFits - - getStandardized - - getVariantIds - - getWeights + - getTraitRuns - title: "Pipelines" desc: > @@ -396,21 +448,21 @@ reference: - autoDecision - raiss - mergeVariantInfo - - mergeSusieCs - filterRelatedness - title: "LD infrastructure" desc: > - Loading and manipulating LD matrices, plus design-matrix - conditioning utilities. + Computing, loading, and manipulating LD matrices, plus + design-matrix conditioning utilities. contents: - loadLdMatrix - loadLdSketch + - loadLdBlock + - computeLd - checkLd - enforceDesignFullRank - ldClumpByScore - ldLoader - - loadLdBlock - ldPruneByCorrelation - filterVariantsByLdReference @@ -421,7 +473,6 @@ reference: lazily, so the Bioconductor accessors apply to it directly. contents: - readGenotypes - - computeLd - loadGenotypeRegion - readAfreq - getRefVariantInfo @@ -467,12 +518,14 @@ reference: - fsusieWrapper - getSusieResult - computeCsTables + - buildTopLoci - extractCsInfo - extractTopPipInfo - formatFinemappingOutput - lbfToAlpha - postprocessFinemappingFits - combineFineMappingResults + - mergeSusieCs - overlapTopLoci - title: "TWAS" @@ -534,11 +587,13 @@ reference: - combineTwasWeights - twasWeightsCv - twasPredict - - twasZ - estimateSparsity - - subtitle: "Multi-method p-value combination" + - subtitle: "Association testing and p-value combination" contents: + # twasZ() turns weights plus GWAS summary statistics into a per-tuple + # Z and p-value, then delegates cross-tuple pooling to combinePValues(). + - twasZ - combinePValues - waldTestPval @@ -653,8 +708,6 @@ reference: - genotypeDelayedArray - as.data.frame.GwasSumStats - rescaleCovW0 - - buildTopLoci - - "getSumStats,GwasSumStats-method" footer: structure: diff --git a/data/gwasSumStatsS4Example.rda b/data/gwasSumStatsS4Example.rda index 18175d97..447dca84 100644 Binary files a/data/gwasSumStatsS4Example.rda and b/data/gwasSumStatsS4Example.rda differ diff --git a/data/multiStudyQtlDatasetExample.rda b/data/multiStudyQtlDatasetExample.rda index 62079697..981d490a 100644 Binary files a/data/multiStudyQtlDatasetExample.rda and b/data/multiStudyQtlDatasetExample.rda differ diff --git a/data/qtlDatasetExample.rda b/data/qtlDatasetExample.rda index a284f4a4..e8cd861a 100644 Binary files a/data/qtlDatasetExample.rda and b/data/qtlDatasetExample.rda differ diff --git a/data/qtlSumStatsExample.rda b/data/qtlSumStatsExample.rda index ef2a0f0c..5f414d64 100644 Binary files a/data/qtlSumStatsExample.rda and b/data/qtlSumStatsExample.rda differ diff --git a/data/qtlSumStatsMulticontextExample.rda b/data/qtlSumStatsMulticontextExample.rda index be6e1c26..fa5b20a8 100644 Binary files a/data/qtlSumStatsMulticontextExample.rda and b/data/qtlSumStatsMulticontextExample.rda differ diff --git a/man/CtwasResult-class.Rd b/man/CtwasResult-class.Rd new file mode 100644 index 00000000..cf833834 --- /dev/null +++ b/man/CtwasResult-class.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/CtwasResult.R +\docType{class} +\name{CtwasResult-class} +\alias{CtwasResult-class} +\title{cTWAS Result Collection} +\description{ +S4 collection of cTWAS runs keyed by the identity tuple + \code{(gwasStudy, study, context, method)}. Each row holds a + \code{\linkS4class{CtwasResultEntry}} payload -- fine-mapping posteriors, + the jointly-estimated group priors, and region metadata -- for one run. +} +\details{ +Unlike the QTL family, \code{trait} is not part of the key: a cTWAS + run is multi-gene, so genes live inside the payload. The optional + \code{jointStudies} / \code{jointContexts} columns tag rows born from a + multi-study or multi-context run and participate in the uniqueness key, + exactly as in the \code{\linkS4class{TwasWeights}} and fine-mapping + families. +} +\seealso{ +\code{\link{CtwasResult}} for the constructor and + \code{\link{ctwasPipeline}} for the pipeline that builds one. +} diff --git a/man/CtwasResultEntry-class.Rd b/man/CtwasResultEntry-class.Rd new file mode 100644 index 00000000..a83a5756 --- /dev/null +++ b/man/CtwasResultEntry-class.Rd @@ -0,0 +1,33 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/CtwasResultEntry.R +\docType{class} +\name{CtwasResultEntry-class} +\alias{CtwasResultEntry-class} +\title{cTWAS Per-Run Payload} +\description{ +Per-run cTWAS payload: the fine-mapping posterior table, the + full per-effect susie alpha table, the jointly-estimated group prior(s), + and per-region metadata. One entry sits in every row of a + \code{\linkS4class{CtwasResult}} collection. +} +\section{Slots}{ + +\describe{ +\item{\code{finemap}}{The per-gene (and, when SNPs are retained, per-SNP) posterior +summary table (\code{ctwas::finemap_regions} \code{finemap_res} shape), +or \code{NULL}.} + +\item{\code{susieAlpha}}{The per-effect susie alpha table +(\code{ctwas::finemap_regions} \code{susie_alpha_res} shape) -- the +fuller cTWAS output retained so the raw run is reconstructable, or +\code{NULL}.} + +\item{\code{param}}{The estimated \code{group_prior} / \code{group_prior_var} for +this run, or \code{NULL}.} + +\item{\code{regionInfo}}{Per-region metadata, or \code{NULL}.} +}} + +\seealso{ +\code{\link{CtwasResultEntry}} for the constructor. +} diff --git a/man/GwasSumStats-class.Rd b/man/GwasSumStats-class.Rd new file mode 100644 index 00000000..c48e31e2 --- /dev/null +++ b/man/GwasSumStats-class.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/gwasSumStats.R +\docType{class} +\name{GwasSumStats-class} +\alias{GwasSumStats-class} +\title{GWAS Summary-Statistic Collection} +\description{ +S4 collection of GWAS summary statistics keyed by the identity + tuple \code{(study)}. Each element is that study's per-variant + \code{GRanges} covering a single LD block, so build one collection per + block when sweeping the genome. +} +\details{ +Required column: \code{study}, unique across rows. The class-level + slots inherited from \code{\linkS4class{SumStatsBase}} -- + \code{ldSketch}, \code{genome} and \code{qcInfo} -- apply uniformly to + every row. +} +\seealso{ +\code{\link{GwasSumStats}} for the constructor and + \code{\linkS4class{QtlSumStats}} for the QTL counterpart. +} diff --git a/man/QtlDataset-class.Rd b/man/QtlDataset-class.Rd index 0d9492eb..d9d45717 100644 --- a/man/QtlDataset-class.Rd +++ b/man/QtlDataset-class.Rd @@ -89,9 +89,7 @@ than one of the things being chosen between, so it is always retained. \describe{ \item{\code{study}}{Character (length 1). Study identifier; used in collection classes to tag downstream \code{FineMappingResult} / \code{TwasWeights} -entries.} - -\item{\code{genotypes}}{The genotype source for lazy access to dosages. +entries. The \code{genotype} experiment's assay reads through this handle; the extraction accessors read it directly, so that QC can be applied per block.} diff --git a/man/QtlSumStats-class.Rd b/man/QtlSumStats-class.Rd new file mode 100644 index 00000000..c909f8f1 --- /dev/null +++ b/man/QtlSumStats-class.Rd @@ -0,0 +1,23 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/qtlSumStats.R +\docType{class} +\name{QtlSumStats-class} +\alias{QtlSumStats-class} +\title{QTL Summary-Statistic Collection} +\description{ +S4 collection of QTL summary statistics keyed by the identity + tuple \code{(study, context, trait)}. Each element is that tuple's + per-variant \code{GRanges} -- \code{variant_id} plus the per-variant + Z / N / MAF mcols. +} +\details{ +Required columns: \code{study}, \code{context} and \code{trait}; + the 3-tuple is unique. The class-level slots inherited from + \code{\linkS4class{SumStatsBase}} -- \code{ldSketch}, \code{genome} and + \code{qcInfo} -- apply uniformly to every row, and \code{qcInfo} records + which \code{\link{summaryStatsQc}} passes have run. +} +\seealso{ +\code{\link{QtlSumStats}} for the constructor and + \code{\linkS4class{GwasSumStats}} for the GWAS counterpart. +} diff --git a/man/RangedTupleList-methods.Rd b/man/RangedTupleList-methods.Rd index bbb59a79..b250d846 100644 --- a/man/RangedTupleList-methods.Rd +++ b/man/RangedTupleList-methods.Rd @@ -10,6 +10,7 @@ \alias{[,RangedTupleList,ANY,ANY,ANY-method} \alias{[[<-,RangedTupleList,ANY,ANY-method} \alias{endoapply,RangedTupleList-method} +\alias{bindROWS,RangedTupleList-method} \title{RangedTupleList methods} \usage{ \S4method{nrow}{RangedTupleList}(x) @@ -27,6 +28,14 @@ \S4method{[[}{RangedTupleList,ANY,ANY}(x, i, j, ...) <- value \S4method{endoapply}{RangedTupleList}(X, FUN, ...) + +\S4method{bindROWS}{RangedTupleList}( + x, + objects = list(), + use.names = TRUE, + ignore.mcols = FALSE, + check = TRUE +) } \arguments{ \item{x, X}{A \code{RangedTupleList} collection.} diff --git a/man/SumStatsBase-class.Rd b/man/SumStatsBase-class.Rd index f70e977a..1f3f2ea5 100644 --- a/man/SumStatsBase-class.Rd +++ b/man/SumStatsBase-class.Rd @@ -8,7 +8,8 @@ Virtual base class for QTL and GWAS summary statistics collections. Concrete subclasses (\code{QtlSumStats}, \code{GwasSumStats}) inherit from \code{\linkS4class{RangedTupleList}} and share the - \code{ldSketch} / \code{genome} / \code{qcInfo} slots. + \code{ldSketch} / \code{qcInfo} slots, and the genome build in + \code{seqinfo()}. Each element is the per-variant \code{GRanges} of one tuple, so \code{x[[i]]} is that tuple's summary statistics and the identity columns @@ -23,8 +24,6 @@ or \code{NULL}. Optional: LD-free workflows (e.g. mash, which operates across conditions per variant) carry \code{NULL}; pipelines that need LD validate its presence when they consume the collection.} -\item{\code{genome}}{Character, genome build label.} - \item{\code{qcInfo}}{A \code{list} recording which QC steps ran. Empty \code{list()} on construction; populated by \code{summaryStatsQc()} with a per-step audit record (filter names, drop counts, liftover target, RAISS settings, etc.). diff --git a/man/TupleRangesView-class.Rd b/man/TupleRangesView-class.Rd deleted file mode 100644 index 4cd731c8..00000000 --- a/man/TupleRangesView-class.Rd +++ /dev/null @@ -1,27 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/TupleRangesView.R -\docType{class} -\name{TupleRangesView-class} -\alias{TupleRangesView-class} -\alias{[,TupleRangesView,ANY,ANY,ANY-method} -\title{Flat view of a tuple collection for plyranges verbs} -\usage{ -\S4method{[}{TupleRangesView,ANY,ANY,ANY}(x, i, j, ..., drop = TRUE) -} -\description{ -A \code{\link[GenomicRanges]{GenomicRanges}} view over a -\code{RangedTupleList}: every element's ranges concatenated, with the -collection's identity tuple broadcast onto each range. Because it is a -\code{GenomicRanges}, plyranges verbs act on it directly. -} -\details{ -Users do not normally build one. The verb methods on \code{RangedTupleList} -convert to a view, delegate, and convert back. -} -\section{Slots}{ - -\describe{ -\item{\code{delegate}}{The flattened \code{GRanges}.} -}} - -\keyword{internal} diff --git a/man/extractCsInfo.Rd b/man/extractCsInfo.Rd index 58a74555..3661559b 100644 --- a/man/extractCsInfo.Rd +++ b/man/extractCsInfo.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/sumstatsQc.R +% Please edit documentation in R/fineMappingWrappers.R \name{extractCsInfo} \alias{extractCsInfo} \title{Process Credible Sets (CS) from Finemapping Results} diff --git a/man/extractTopPipInfo.Rd b/man/extractTopPipInfo.Rd index 2c1c3b77..6851486e 100644 --- a/man/extractTopPipInfo.Rd +++ b/man/extractTopPipInfo.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/sumstatsQc.R +% Please edit documentation in R/fineMappingWrappers.R \name{extractTopPipInfo} \alias{extractTopPipInfo} \title{Extract Information for Top Variant from Finemapping Results} diff --git a/man/flattenTupleRanges.Rd b/man/flattenTupleRanges.Rd index ff4282be..f2c1f81c 100644 --- a/man/flattenTupleRanges.Rd +++ b/man/flattenTupleRanges.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/TupleRangesView.R +% Please edit documentation in R/RangedTupleList.R \name{flattenTupleRanges} \alias{flattenTupleRanges} \title{Flatten a tuple collection to a plain GRanges} diff --git a/man/getSusieResult.Rd b/man/getSusieResult.Rd index 31cc76f7..928f5060 100644 --- a/man/getSusieResult.Rd +++ b/man/getSusieResult.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/sumstatsQc.R +% Please edit documentation in R/fineMappingWrappers.R \name{getSusieResult} \alias{getSusieResult} \title{Extract the trimmed SuSiE fit from a finemapping pipeline result} diff --git a/man/nestTupleRanges.Rd b/man/nestTupleRanges.Rd index 95ddf771..54d19155 100644 --- a/man/nestTupleRanges.Rd +++ b/man/nestTupleRanges.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/TupleRangesView.R +% Please edit documentation in R/RangedTupleList.R \name{nestTupleRanges} \alias{nestTupleRanges} \title{Re-nest a flattened GRanges into a tuple collection} diff --git a/man/show-methods.Rd b/man/show-methods.Rd index 0c54b6fb..63613ad7 100644 --- a/man/show-methods.Rd +++ b/man/show-methods.Rd @@ -3,9 +3,8 @@ % R/AnnotationMatrix.R, R/ColocResult.R, R/ColocBoostResult.R, % R/CtwasResult.R, R/GwasFineMappingResult.R, R/H2Estimate.R, R/LdData.R, % R/LdEigen.R, R/LdScore.R, R/MashPrior.R, R/QtlDataset.R, -% R/MultiStudyQtlDataset.R, R/QtlFineMappingResult.R, R/TupleRangesView.R, -% R/qtlSumStats.R, R/gwasSumStats.R, R/fineMappingRow.R, R/twasWeights.R, -% R/twasWeightsRow.R +% R/MultiStudyQtlDataset.R, R/QtlFineMappingResult.R, R/qtlSumStats.R, +% R/gwasSumStats.R, R/fineMappingRow.R, R/twasWeights.R, R/twasWeightsRow.R \name{show-methods} \alias{show-methods} \alias{show,GenotypeHandle-method} @@ -22,7 +21,6 @@ \alias{show,QtlDataset-method} \alias{show,MultiStudyQtlDataset-method} \alias{show,QtlFineMappingResult-method} -\alias{show,TupleRangesView-method} \alias{show,QtlSumStats-method} \alias{show,GwasSumStats-method} \alias{show,FineMappingRow-method} @@ -58,8 +56,6 @@ \S4method{show}{QtlFineMappingResult}(object) -\S4method{show}{TupleRangesView}(object) - \S4method{show}{QtlSumStats}(object) \S4method{show}{GwasSumStats}(object) diff --git a/tests/testthat/test_RangedTupleList.R b/tests/testthat/test_RangedTupleList.R index 64cce75f..2b3d9e18 100644 --- a/tests/testthat/test_RangedTupleList.R +++ b/tests/testthat/test_RangedTupleList.R @@ -562,3 +562,255 @@ test_that(".rtlGatherElements does not re-validate the collection", { x <- .rtl_makeKid() expect_identical(.rtlGatherElements(x, 2L), x[[2]]) }) + + +# =========================================================================== +# The plyranges bridge: the flatten / nest conversions and the dplyr verbs +# built on them. +# +# plyranges has no GRangesList support, so the collections cannot be operated +# on directly. The bridge flattens to a plain GenomicRanges -- which plyranges +# does support natively -- and nests the result back again. +# =========================================================================== + +test_that("flattenTupleRanges concatenates every element's ranges", { + data(qtlSumStatsMulticontextExample, envir = environment()) + mc <- qtlSumStatsMulticontextExample + flat <- flattenTupleRanges(mc) + expect_s4_class(flat, "GRanges") + expect_equal(length(flat), sum(lengths(mc))) +}) + +test_that("flattenTupleRanges broadcasts the identity tuple, dot-prefixed", { + data(qtlSumStatsMulticontextExample, envir = environment()) + mc <- qtlSumStatsMulticontextExample + flat <- flattenTupleRanges(mc) + cn <- colnames(mcols(flat)) + expect_true(all(is_in(c(".study", ".context", ".trait"), cn))) + # One value per range, matching the element it came from. + expect_setequal(unique(mcols(flat)$.context), as.character(mc$context)) + expect_equal( + sum(mcols(flat)$.context == "blood"), + lengths(mc)[[which(mc$context == "blood")]] + ) +}) + +test_that("the dot prefix keeps a colliding per-range column intact", { + # On gwasFineMappingExample the outer `method` is "susie" while the + # per-range `method` is "susieRss" -- the collection's label against the + # fitter actually used. Broadcasting onto the bare name would overwrite + # one with the other. + data(gwasFineMappingExample, envir = environment()) + x <- gwasFineMappingExample + expect_true(is_in("method", colnames(mcols(x[[1]])))) + flat <- flattenTupleRanges(x) + expect_true(all(is_in(c("method", ".method"), colnames(mcols(flat))))) + expect_false(identical( + unique(mcols(flat)$method), + unique(mcols(flat)$.method) + )) +}) + +test_that("flattenTupleRanges drops non-broadcastable payload columns", { + # susieFit / cvResult describe an ELEMENT, so there is no range to put + # them on. + data(qtlFineMappingExample, envir = environment()) + flat <- flattenTupleRanges(qtlFineMappingExample) + expect_false(is_in("susieFit", colnames(mcols(flat)))) + expect_false(is_in(".susieFit", colnames(mcols(flat)))) + expect_true(is_in(".study", colnames(mcols(flat)))) +}) + +test_that("flattenTupleRanges rejects a non-collection", { + expect_error( + flattenTupleRanges(GenomicRanges::GRanges()), + "RangedTupleList" + ) +}) + +test_that("nestTupleRanges round-trips an unmodified flatten", { + data(qtlSumStatsMulticontextExample, envir = environment()) + mc <- qtlSumStatsMulticontextExample + back <- nestTupleRanges(flattenTupleRanges(mc), mc) + expect_s4_class(back, "QtlSumStats") + expect_equal(lengths(back), lengths(mc)) + expect_equal(mcols(back), mcols(mc)) + # Collection-level slots survive, which is the point of nesting through a + # template rather than rebuilding from the ranges alone. + expect_identical(getLdSketch(back), getLdSketch(mc)) + expect_identical(getGenome(back), getGenome(mc)) +}) + +test_that("nestTupleRanges returns a filtered flatten to its elements", { + data(qtlSumStatsMulticontextExample, envir = environment()) + mc <- qtlSumStatsMulticontextExample + flat <- flattenTupleRanges(mc) + kept <- flat[mcols(flat)$Z > 0] + back <- nestTupleRanges(kept, mc) + expect_equal(sum(lengths(back)), length(kept)) + expect_lt(sum(lengths(back)), sum(lengths(mc))) + # Every range lands under the tuple it carried. + idx <- which(mc$context == "blood") + expect_equal( + lengths(back)[[idx]], + sum(mcols(kept)$.context == "blood") + ) +}) + +test_that("nestTupleRanges keeps a tuple whose ranges all vanished", { + # An empty element, not a dropped row: the collection keeps its shape so + # its metadata stays aligned. + data(qtlSumStatsMulticontextExample, envir = environment()) + mc <- qtlSumStatsMulticontextExample + flat <- flattenTupleRanges(mc) + onlyBlood <- flat[mcols(flat)$.context == "blood"] + back <- nestTupleRanges(onlyBlood, mc) + expect_equal(nrow(back), nrow(mc)) + expect_equal(sum(lengths(back) == 0L), nrow(mc) - 1L) +}) + +test_that("nestTupleRanges survives an NA-valued identity column", { + # str_c() propagates NA, so an NA identity value (varY is NA_real_ on a + # z-score collection) would make every key NA and the row lookup a + # logical subscript full of NAs. + data(qtlSumStatsMulticontextExample, envir = environment()) + mc <- qtlSumStatsMulticontextExample + expect_true(any(is.na(mcols(mc)$varY))) + expect_equal( + lengths(nestTupleRanges(flattenTupleRanges(mc), mc)), + lengths(mc) + ) +}) + +test_that("nestTupleRanges validates its arguments", { + data(qtlSumStatsMulticontextExample, envir = environment()) + mc <- qtlSumStatsMulticontextExample + flat <- flattenTupleRanges(mc) + expect_error(nestTupleRanges(flat, "not a collection"), "RangedTupleList") + expect_error(nestTupleRanges("not ranges", mc), "GRanges") +}) + + + + +# --------------------------------------------------------------------------- +# plyranges / dplyr verbs. +# +# These call the methods directly rather than through the generic: dispatch +# needs the S3method() entries a roxygen regen writes, and the logic under test +# is the same either way. The final test checks dispatch itself, and skips +# until the registration exists. +# --------------------------------------------------------------------------- + +# @noRd +.trv_mc <- function() { + data(qtlSumStatsMulticontextExample, envir = environment()) + qtlSumStatsMulticontextExample +} + +test_that("filter keeps the collection and narrows its elements", { + mc <- .trv_mc() + out <- filter.RangedTupleList(mc, Z > 2) + expect_s4_class(out, "QtlSumStats") + expect_equal(nrow(out), nrow(mc)) + expect_lt(sum(lengths(out)), sum(lengths(mc))) + expect_identical(getLdSketch(out), getLdSketch(mc)) +}) + +test_that("filter can mix identity and per-range columns", { + # The whole point of broadcasting the tuple: one predicate spanning both. + mc <- .trv_mc() + out <- filter.RangedTupleList(mc, .context == "blood" & Z > 1) + idx <- which(mc$context == "blood") + expect_gt(lengths(out)[[idx]], 0L) + expect_equal(sum(lengths(out)[-idx]), 0L) +}) + +test_that("mutate adds a per-range column without changing the shape", { + mc <- .trv_mc() + out <- mutate.RangedTupleList(mc, hit = Z > 2) + expect_equal(lengths(out), lengths(mc)) + expect_true(is_in("hit", colnames(mcols(out[[1]])))) +}) + +test_that("select keeps the identity columns it needs to nest", { + # Dropping them would leave the ranges with no tuple, and every one would + # fall into the first element. + mc <- .trv_mc() + out <- select.RangedTupleList(mc, Z) + expect_equal(lengths(out), lengths(mc)) +}) + +test_that("slice takes n ranges from EACH element", { + # Slicing the flattened set would take n in total and empty all but the + # first element. + mc <- .trv_mc() + out <- slice.RangedTupleList(mc, 1:5) + expect_equal(unname(lengths(out)), rep(5L, nrow(mc))) +}) + +test_that("arrange orders within each element", { + mc <- .trv_mc() + out <- arrange.RangedTupleList(mc, Z) + expect_equal(lengths(out), lengths(mc)) + for (i in seq_len(nrow(out))) { + expect_false(is.unsorted(mcols(out[[i]])$Z)) + } +}) + +test_that("summarise reduces to a table rather than a collection", { + mc <- .trv_mc() + out <- summarise.RangedTupleList(mc, n = plyranges::n()) + expect_false(is(out, "RangedTupleList")) + expect_equal(NROW(out), 1L) +}) + +test_that("group_by returns plyranges' own grouped representation", { + mc <- .trv_mc() + out <- group_by.RangedTupleList(mc, .context) + expect_s4_class(out, "GroupedGenomicRanges") + expect_equal(length(out), sum(lengths(mc))) +}) + + +test_that("the verbs are reachable through the dplyr generics", { + # Needs the S3method() entries from a roxygen regen; skips until then. + # Probed by attempting dispatch rather than by looking the function up: + # getS3method() finds it in the namespace whether or not it is registered. + mc <- .trv_mc() + dispatches <- tryCatch( + { + dplyr::filter(mc, Z > 2) + TRUE + }, + error = function(e) FALSE + ) + skip_if_not(dispatches, "S3 methods not yet registered (regen pending)") + expect_s4_class(dplyr::filter(mc, Z > 2), "QtlSumStats") + expect_s4_class(dplyr::slice(mc, 1:2), "QtlSumStats") +}) + + +# =========================================================================== +# Degenerate collections and the plyranges requirement +# =========================================================================== + +test_that("flattening a collection with no identity columns keeps its ranges", { + # With nothing to key on, every range belongs to the single element -- + # which is what a one-row collection means. + mc <- .trv_mc() + bare <- mc[1, ] + mcols(bare) <- mcols(bare)[, character(0), drop = FALSE] + flat <- flattenTupleRanges(bare) + expect_equal(length(flat), sum(lengths(bare))) +}) + +test_that(".rtlTupleKeys gives every row the same key with no key columns", { + md <- S4Vectors::DataFrame(entry = S4Vectors::SimpleList(1, 2)) + expect_equal(.rtlTupleKeys(md, character(0), 3L), rep("", 3L)) +}) + +test_that(".rtlTupleKeyCols is empty when there are no mcols", { + gr <- GenomicRanges::GRanges("chr1", IRanges::IRanges(1, 2)) + expect_equal(.rtlTupleKeyCols(gr), character(0)) +}) diff --git a/tests/testthat/test_TupleRangesView.R b/tests/testthat/test_TupleRangesView.R deleted file mode 100644 index 58c3850c..00000000 --- a/tests/testthat/test_TupleRangesView.R +++ /dev/null @@ -1,314 +0,0 @@ -# Tests for the plyranges bridge: TupleRangesView and the flatten / nest -# conversions. -# -# plyranges has no GRangesList support, so the collections cannot be operated -# on directly. The bridge converts to a flat GenomicRanges view -- which -# plyranges does support, via the DelegatingGenomicRanges extension point -- -# and back again. - -test_that("flattenTupleRanges concatenates every element's ranges", { - data(qtlSumStatsMulticontextExample, envir = environment()) - mc <- qtlSumStatsMulticontextExample - flat <- flattenTupleRanges(mc) - expect_s4_class(flat, "GRanges") - expect_equal(length(flat), sum(lengths(mc))) -}) - -test_that("flattenTupleRanges broadcasts the identity tuple, dot-prefixed", { - data(qtlSumStatsMulticontextExample, envir = environment()) - mc <- qtlSumStatsMulticontextExample - flat <- flattenTupleRanges(mc) - cn <- colnames(mcols(flat)) - expect_true(all(is_in(c(".study", ".context", ".trait"), cn))) - # One value per range, matching the element it came from. - expect_setequal(unique(mcols(flat)$.context), as.character(mc$context)) - expect_equal( - sum(mcols(flat)$.context == "blood"), - lengths(mc)[[which(mc$context == "blood")]] - ) -}) - -test_that("the dot prefix keeps a colliding per-range column intact", { - # On gwasFineMappingExample the outer `method` is "susie" while the - # per-range `method` is "susieRss" -- the collection's label against the - # fitter actually used. Broadcasting onto the bare name would overwrite - # one with the other. - data(gwasFineMappingExample, envir = environment()) - x <- gwasFineMappingExample - expect_true(is_in("method", colnames(mcols(x[[1]])))) - flat <- flattenTupleRanges(x) - expect_true(all(is_in(c("method", ".method"), colnames(mcols(flat))))) - expect_false(identical( - unique(mcols(flat)$method), - unique(mcols(flat)$.method) - )) -}) - -test_that("flattenTupleRanges drops non-broadcastable payload columns", { - # susieFit / cvResult describe an ELEMENT, so there is no range to put - # them on. - data(qtlFineMappingExample, envir = environment()) - flat <- flattenTupleRanges(qtlFineMappingExample) - expect_false(is_in("susieFit", colnames(mcols(flat)))) - expect_false(is_in(".susieFit", colnames(mcols(flat)))) - expect_true(is_in(".study", colnames(mcols(flat)))) -}) - -test_that("flattenTupleRanges rejects a non-collection", { - expect_error( - flattenTupleRanges(GenomicRanges::GRanges()), - "RangedTupleList" - ) -}) - -test_that("nestTupleRanges round-trips an unmodified flatten", { - data(qtlSumStatsMulticontextExample, envir = environment()) - mc <- qtlSumStatsMulticontextExample - back <- nestTupleRanges(flattenTupleRanges(mc), mc) - expect_s4_class(back, "QtlSumStats") - expect_equal(lengths(back), lengths(mc)) - expect_equal(mcols(back), mcols(mc)) - # Collection-level slots survive, which is the point of nesting through a - # template rather than rebuilding from the ranges alone. - expect_identical(getLdSketch(back), getLdSketch(mc)) - expect_identical(getGenome(back), getGenome(mc)) -}) - -test_that("nestTupleRanges returns a filtered flatten to its elements", { - data(qtlSumStatsMulticontextExample, envir = environment()) - mc <- qtlSumStatsMulticontextExample - flat <- flattenTupleRanges(mc) - kept <- flat[mcols(flat)$Z > 0] - back <- nestTupleRanges(kept, mc) - expect_equal(sum(lengths(back)), length(kept)) - expect_lt(sum(lengths(back)), sum(lengths(mc))) - # Every range lands under the tuple it carried. - idx <- which(mc$context == "blood") - expect_equal( - lengths(back)[[idx]], - sum(mcols(kept)$.context == "blood") - ) -}) - -test_that("nestTupleRanges keeps a tuple whose ranges all vanished", { - # An empty element, not a dropped row: the collection keeps its shape so - # its metadata stays aligned. - data(qtlSumStatsMulticontextExample, envir = environment()) - mc <- qtlSumStatsMulticontextExample - flat <- flattenTupleRanges(mc) - onlyBlood <- flat[mcols(flat)$.context == "blood"] - back <- nestTupleRanges(onlyBlood, mc) - expect_equal(nrow(back), nrow(mc)) - expect_equal(sum(lengths(back) == 0L), nrow(mc) - 1L) -}) - -test_that("nestTupleRanges survives an NA-valued identity column", { - # str_c() propagates NA, so an NA identity value (varY is NA_real_ on a - # z-score collection) would make every key NA and the row lookup a - # logical subscript full of NAs. - data(qtlSumStatsMulticontextExample, envir = environment()) - mc <- qtlSumStatsMulticontextExample - expect_true(any(is.na(mcols(mc)$varY))) - expect_equal( - lengths(nestTupleRanges(flattenTupleRanges(mc), mc)), - lengths(mc) - ) -}) - -test_that("nestTupleRanges validates its arguments", { - data(qtlSumStatsMulticontextExample, envir = environment()) - mc <- qtlSumStatsMulticontextExample - flat <- flattenTupleRanges(mc) - expect_error(nestTupleRanges(flat, "not a collection"), "RangedTupleList") - expect_error(nestTupleRanges("not ranges", mc), "GRanges") -}) - -test_that("a TupleRangesView is a GenomicRanges and subsets to itself", { - data(qtlSumStatsMulticontextExample, envir = environment()) - flat <- flattenTupleRanges(qtlSumStatsMulticontextExample) - v <- .tupleRangesView(flat) - expect_s4_class(v, "TupleRangesView") - # Being a GenomicRanges is what makes plyranges dispatch at all. - expect_true(is(v, "GenomicRanges")) - expect_equal(length(v), length(flat)) - expect_s4_class(v[1:10], "TupleRangesView") - expect_equal(length(v[1:10]), 10L) -}) - -test_that("the view constructor keeps mcols parallel to the object", { - # DelegatingGenomicRanges carries its own elementMetadata; leaving it - # zero-row makes the class reject itself with "'mcols(x)' is not parallel - # to 'x'". - data(qtlSumStatsMulticontextExample, envir = environment()) - flat <- flattenTupleRanges(qtlSumStatsMulticontextExample) - expect_silent(v <- .tupleRangesView(flat)) - expect_true(validObject(v)) -}) - -test_that("the view rejects a two-argument subscript", { - data(qtlSumStatsMulticontextExample, envir = environment()) - v <- .tupleRangesView(flattenTupleRanges(qtlSumStatsMulticontextExample)) - expect_error(v[1, 1], "two-argument") -}) - - -# --------------------------------------------------------------------------- -# plyranges / dplyr verbs. -# -# These call the methods directly rather than through the generic: dispatch -# needs the S3method() entries a roxygen regen writes, and the logic under test -# is the same either way. The final test checks dispatch itself, and skips -# until the registration exists. -# --------------------------------------------------------------------------- - -# @noRd -.trv_mc <- function() { - data(qtlSumStatsMulticontextExample, envir = environment()) - qtlSumStatsMulticontextExample -} - -test_that("filter keeps the collection and narrows its elements", { - mc <- .trv_mc() - out <- filter.RangedTupleList(mc, Z > 2) - expect_s4_class(out, "QtlSumStats") - expect_equal(nrow(out), nrow(mc)) - expect_lt(sum(lengths(out)), sum(lengths(mc))) - expect_identical(getLdSketch(out), getLdSketch(mc)) -}) - -test_that("filter can mix identity and per-range columns", { - # The whole point of broadcasting the tuple: one predicate spanning both. - mc <- .trv_mc() - out <- filter.RangedTupleList(mc, .context == "blood" & Z > 1) - idx <- which(mc$context == "blood") - expect_gt(lengths(out)[[idx]], 0L) - expect_equal(sum(lengths(out)[-idx]), 0L) -}) - -test_that("mutate adds a per-range column without changing the shape", { - mc <- .trv_mc() - out <- mutate.RangedTupleList(mc, hit = Z > 2) - expect_equal(lengths(out), lengths(mc)) - expect_true(is_in("hit", colnames(mcols(out[[1]])))) -}) - -test_that("select keeps the identity columns it needs to nest", { - # Dropping them would leave the ranges with no tuple, and every one would - # fall into the first element. - mc <- .trv_mc() - out <- select.RangedTupleList(mc, Z) - expect_equal(lengths(out), lengths(mc)) -}) - -test_that("slice takes n ranges from EACH element", { - # Slicing the flattened set would take n in total and empty all but the - # first element. - mc <- .trv_mc() - out <- slice.RangedTupleList(mc, 1:5) - expect_equal(unname(lengths(out)), rep(5L, nrow(mc))) -}) - -test_that("arrange orders within each element", { - mc <- .trv_mc() - out <- arrange.RangedTupleList(mc, Z) - expect_equal(lengths(out), lengths(mc)) - for (i in seq_len(nrow(out))) { - expect_false(is.unsorted(mcols(out[[i]])$Z)) - } -}) - -test_that("summarise reduces to a table rather than a collection", { - mc <- .trv_mc() - out <- summarise.RangedTupleList(mc, n = plyranges::n()) - expect_false(is(out, "RangedTupleList")) - expect_equal(NROW(out), 1L) -}) - -test_that("group_by returns plyranges' own grouped representation", { - mc <- .trv_mc() - out <- group_by.RangedTupleList(mc, .context) - expect_s4_class(out, "GroupedGenomicRanges") - expect_equal(length(out), sum(lengths(mc))) -}) - -test_that("select and group_by work on the view itself", { - # plyranges' own DelegatingGenomicRanges support misses both: they work on - # a plain GRanges and fail on any delegating class. These patch that. - mc <- .trv_mc() - v <- .tupleRangesView(flattenTupleRanges(mc)) - expect_s4_class(select.TupleRangesView(v, Z), "TupleRangesView") - expect_s4_class( - group_by.TupleRangesView(v, .context), - "GroupedGenomicRanges" - ) -}) - -test_that("the verbs are reachable through the dplyr generics", { - # Needs the S3method() entries from a roxygen regen; skips until then. - # Probed by attempting dispatch rather than by looking the function up: - # getS3method() finds it in the namespace whether or not it is registered. - mc <- .trv_mc() - dispatches <- tryCatch( - { - dplyr::filter(mc, Z > 2) - TRUE - }, - error = function(e) FALSE - ) - skip_if_not(dispatches, "S3 methods not yet registered (regen pending)") - expect_s4_class(dplyr::filter(mc, Z > 2), "QtlSumStats") - expect_s4_class(dplyr::slice(mc, 1:2), "QtlSumStats") -}) - -test_that("show reports the view's range and metadata-column counts", { - # The view is internal, but its show method is what a developer sees when - # one surfaces in a browser() or an error trace. - data(qtlSumStatsMulticontextExample, envir = environment()) - flat <- flattenTupleRanges(qtlSumStatsMulticontextExample) - v <- .tupleRangesView(flat) - expect_output(show(v), "TupleRangesView:") - expect_output(show(v), str_c(length(flat), " ranges")) - expect_output(show(v), "metadata column\\(s\\)") -}) - -# =========================================================================== -# Degenerate collections and the plyranges requirement -# =========================================================================== - -test_that("flattening a collection with no identity columns keeps its ranges", { - # With nothing to key on, every range belongs to the single element -- - # which is what a one-row collection means. - mc <- .trv_mc() - bare <- mc[1, ] - mcols(bare) <- mcols(bare)[, character(0), drop = FALSE] - flat <- flattenTupleRanges(bare) - expect_equal(length(flat), sum(lengths(bare))) -}) - -test_that(".rtlTupleKeys gives every row the same key with no key columns", { - md <- S4Vectors::DataFrame(entry = S4Vectors::SimpleList(1, 2)) - expect_equal(.rtlTupleKeys(md, character(0), 3L), rep("", 3L)) -}) - -test_that(".rtlTupleKeyCols is empty when there are no mcols", { - gr <- GenomicRanges::GRanges("chr1", IRanges::IRanges(1, 2)) - expect_equal(.rtlTupleKeyCols(gr), character(0)) -}) - -test_that(".rtlAsRanges unwraps a view and passes a GRanges through", { - flat <- flattenTupleRanges(.trv_mc()) - v <- .tupleRangesView(flat) - expect_s4_class(.rtlAsRanges(v), "GRanges") - expect_false(methods::is(.rtlAsRanges(v), "TupleRangesView")) - expect_identical(.rtlAsRanges(flat), flat) -}) - -test_that("select(.drop_ranges = TRUE) returns the bare table, not a view", { - # plyranges' own escape hatch: the caller asked for a data frame rather - # than ranges, so re-wrapping it as a view would undo the request. - flat <- flattenTupleRanges(.trv_mc()) - v <- .tupleRangesView(flat) - out <- select.TupleRangesView(v, Z, .drop_ranges = TRUE) - expect_false(methods::is(out, "TupleRangesView")) - expect_false(methods::is(out, "GRanges")) -}) diff --git a/tests/testthat/test_ctwasPipeline.R b/tests/testthat/test_ctwasPipeline.R index e7d4ebff..27daa80b 100644 --- a/tests/testthat/test_ctwasPipeline.R +++ b/tests/testthat/test_ctwasPipeline.R @@ -1642,7 +1642,7 @@ test_that("ctwasPipeline: dispatches assemble → est → screen → finemap and expect_equal(as.character(out$gwasStudy), "G1") # The run's jointly-estimated param is carried on the row. expect_equal( - unname(getCtwasParam(pecotmr:::.collectionEntry(out, 1L))$group_prior), + unname(getCtwasParam(out$entry[[1L]])$group_prior), c(0.1, 0.0001) ) # getFinemap aggregates the per-gene rows, tagged with run identity. @@ -2147,7 +2147,7 @@ test_that("ctwasPipeline: real-engine end-to-end on the bundled example panel", data(qtlDatasetExample) gss <- gwasSumStatsS4Example qd <- qtlDatasetExample - gh <- qd@genotypes + gh <- getGenotypeHandle(qd) # Two 5-variant synthetic genes from the bundled panel, one anchored in # each LD block below. cTWAS's EM runs per region and cannot fit a block @@ -2353,7 +2353,7 @@ test_that(".ctwasFilterVariants: returns NULL when no variants survive", { test_that(".ctwasBuildWeights: maxNumVariants caps the per-gene weight matrix", { data(qtlDatasetExample) qd <- qtlDatasetExample - gh <- qd@genotypes + gh <- getGenotypeHandle(qd) vids <- getSnpInfo(gh)$SNP[1:5] ent <- twasWeightsRow( variantIds = vids, @@ -2378,7 +2378,7 @@ test_that(".ctwasBuildWeights: maxNumVariants caps the per-gene weight matrix", test_that(".ctwasBuildWeights: twasWeightCutoff drops low-magnitude variants", { data(qtlDatasetExample) qd <- qtlDatasetExample - gh <- qd@genotypes + gh <- getGenotypeHandle(qd) vids <- getSnpInfo(gh)$SNP[1:5] ent <- twasWeightsRow( variantIds = vids, diff --git a/tests/testthat/test_fineMappingWrappers.R b/tests/testthat/test_fineMappingWrappers.R index 7131a23d..33be971e 100644 --- a/tests/testthat/test_fineMappingWrappers.R +++ b/tests/testthat/test_fineMappingWrappers.R @@ -3216,3 +3216,225 @@ test_that(".fmEffectIndices parses L-names and falls back to position", { expect_equal(.fmEffectIndices(list(L1 = 1, L2 = 2, L7 = 3)), c(1L, 2L, 7L)) expect_equal(.fmEffectIndices(list(1, 2, 3)), c(1L, 2L, 3L)) }) + + +context("post_finemapping_cs_extraction") + +# A single-row QtlFineMappingResult standing in for the retired +# FineMappingRow. topLoci defaults to one row per variant: the row-payload +# builder requires the two to be aligned, where the entry tolerated an empty +# table beside a non-empty variant list. +.testFineMappingRow <- function( + variantIds, + susieFit = list(), + topLoci = NULL +) { + if (is.null(topLoci)) { + topLoci <- data.frame( + variant_id = variantIds, + pip = rep(0, length(variantIds)), + stringsAsFactors = FALSE + ) + } + QtlFineMappingResult( + study = "s1", + context = "c1", + trait = "t1", + method = "susie", + entry = list(fineMappingRow( + variantIds = variantIds, + susieFit = susieFit, + topLoci = topLoci + )) + ) +} + +# =========================================================================== +# getSusieResult +# =========================================================================== + +test_that("getSusieResult returns NULL for empty input", { + result <- getSusieResult(list()) + expect_null(result) +}) + +test_that("getSusieResult returns NULL when finemappingEntry missing", { + result <- getSusieResult(list(some_data = 42)) + expect_null(result) +}) + +test_that("getSusieResult returns trimmed result when present", { + mock_result <- list(pip = c(0.1, 0.5, 0.3), sets = list(cs = list())) + con_data <- list( + finemappingEntry = .testFineMappingRow( + variantIds = c("chr1:100:A:G", "chr1:200:C:T", "chr1:300:G:A"), + susieFit = mock_result + ) + ) + result <- getSusieResult(con_data) + expect_equal(result, mock_result) +}) + +# =========================================================================== +# extractTopPipInfo +# =========================================================================== + +test_that("extractTopPipInfo finds top PIP variant", { + con_data <- list( + finemappingEntry = .testFineMappingRow( + variantIds = c("chr1:100:A:G", "chr1:200:C:T", "chr1:300:G:A"), + susieFit = list(pip = c(0.1, 0.7, 0.2)) + ), + sumstats = list(z = c(1.0, 3.5, -0.5)) + ) + result <- extractTopPipInfo(con_data$finemappingEntry, con_data$sumstats) + expect_equal(result$top_variant, "chr1:200:C:T") + expect_equal(result$top_pip, 0.7) + expect_equal(result$top_z, 3.5) + expect_equal(result$top_variant_index, 2) + expect_true(is.na(result$cs_name)) + expect_true(is.na(result$variants_per_cs)) +}) + +test_that("extractTopPipInfo computes p_value from z", { + con_data <- list( + finemappingEntry = .testFineMappingRow( + variantIds = c("chr1:100:A:G", "chr1:200:C:T", "chr1:300:G:A"), + susieFit = list(pip = c(0.9, 0.05, 0.05)) + ), + sumstats = list(z = c(5.0, 0.5, -0.3)) + ) + result <- extractTopPipInfo(con_data$finemappingEntry, con_data$sumstats) + expected_pval <- pecotmr:::.zToPvalue(5.0) + expect_equal(result$p_value, expected_pval) +}) + +test_that("extractTopPipInfo handles ties by taking first max", { + con_data <- list( + finemappingEntry = .testFineMappingRow( + variantIds = c("chr1:100:A:G", "chr1:200:C:T", "chr1:300:G:A"), + susieFit = list(pip = c(0.5, 0.5, 0.5)) + ), + sumstats = list(z = c(1.0, 2.0, 3.0)) + ) + result <- extractTopPipInfo(con_data$finemappingEntry, con_data$sumstats) + expect_equal(result$top_variant_index, 1) + expect_equal(result$top_pip, 0.5) +}) + +# =========================================================================== +# extractCsInfo +# =========================================================================== + +test_that("extractCsInfo extracts single CS correctly", { + data(qtlSumStatsExample) + fe <- .testFineMappingRow( + variantIds = c("chr1:100:A:G", "chr1:200:C:T", "chr1:300:G:A"), + susieFit = list(sets = list(cs = list(L_1 = c(1, 2)))) + ) + top_loci_table <- data.frame( + variant_id = c("chr1:100:A:G", "chr1:200:C:T"), + pip = c(0.3, 0.8), + z = c(2.0, 4.5), + stringsAsFactors = FALSE + ) + # A single CS short-circuits (no between-CS correlation), so the unrelated + # ldSource is not consulted. + result <- extractCsInfo( + fe, + csNames = "L_1", + topLociTable = top_loci_table, + ldSource = qtlSumStatsExample + ) + expect_equal(nrow(result), 1) + expect_equal(result$cs_name, "L_1") + expect_equal(result$top_variant, "chr1:200:C:T") + expect_equal(result$top_pip, 0.8) + expect_equal(result$variants_per_cs, 2) + expect_true(is.na(result$cs_corr_max)) + expect_true(is.na(result$cs_corr_min)) + expect_false("cs_corr_1" %in% colnames(result)) +}) + +test_that("extractCsInfo builds correlation columns from the ldSource", { + data(qtlSumStatsExample) + ss <- qtlSumStatsExample + vids <- rownames(getLdSketch(ss)) + set.seed(1) + fit <- list( + sets = list(cs = list(L_1 = c(1L, 2L, 3L), L_2 = c(90L, 91L))), + pip = runif(length(vids)) + ) + tl <- data.frame( + variant_id = vids, + pip = fit$pip, + z = rnorm(length(vids)), + stringsAsFactors = FALSE + ) + fe <- .testFineMappingRow( + variantIds = vids, + susieFit = fit, + topLoci = tl + ) + result <- extractCsInfo( + fe, + csNames = c("L_1", "L_2"), + topLociTable = tl, + ldSource = ss + ) + expect_equal(nrow(result), 2) + expect_true(all( + c("cs_corr_1", "cs_corr_2", "cs_corr_max", "cs_corr_min") %in% + colnames(result) + )) + # The cs_corr_j columns are the columns of the computed between-CS matrix + # (symmetric; diagonal == 1), reduced on demand from the ldSource. + cc <- computeCsCorrelation(fe, ss) + expect_equal(result$cs_corr_1, unname(cc[, 1])) + expect_equal(result$cs_corr_2, unname(cc[, 2])) + expect_equal(result$cs_corr_max, rep(abs(cc[1, 2]), 2)) + expect_equal(result$cs_corr_min, rep(abs(cc[1, 2]), 2)) +}) + +test_that("extractCsInfo computes p_value from z-score", { + data(qtlSumStatsExample) + fe <- .testFineMappingRow( + variantIds = c("chr1:100:A:G", "chr1:200:C:T"), + susieFit = list(sets = list(cs = list(L_1 = c(1, 2)))) + ) + top_loci_table <- data.frame( + variant_id = c("chr1:100:A:G", "chr1:200:C:T"), + pip = c(0.9, 0.1), + z = c(5.0, 0.5), + stringsAsFactors = FALSE + ) + result <- extractCsInfo( + fe, + csNames = "L_1", + topLociTable = top_loci_table, + ldSource = qtlSumStatsExample + ) + expected_pval <- pecotmr:::.zToPvalue(5.0) + expect_equal(result$p_value, expected_pval, tolerance = 1e-10) +}) + +# =========================================================================== +# getSusieResult: trimmed fit is empty +# =========================================================================== + +test_that("getSusieResult returns NULL when the trimmed susie fit is empty", { + conData <- list( + # An empty fit with variants present: topLoci must still be aligned + # row-for-row, which the entry did not enforce. + finemappingEntry = fineMappingRow( + variantIds = c("chr1:100:A:G", "chr1:200:C:T"), + susieFit = list(), + topLoci = data.frame( + variant_id = c("chr1:100:A:G", "chr1:200:C:T"), + pip = c(0, 0), + stringsAsFactors = FALSE + ) + ) + ) + expect_null(getSusieResult(conData)) +}) diff --git a/tests/testthat/test_ld.R b/tests/testthat/test_ld.R index cf059f71..376aeb48 100644 --- a/tests/testthat/test_ld.R +++ b/tests/testthat/test_ld.R @@ -4463,14 +4463,14 @@ test_that(".panelVariantStats reports NA MAF for an all-missing variant", { test_that(".panelVariantFilter is a no-op at its defaults", { data(qtlDatasetExample) gh <- getGenotypes(qtlDatasetExample) - handle <- qtlDatasetExample@genotypes + handle <- getGenotypeHandle(qtlDatasetExample) ids <- normalizeVariantId(getSnpInfo(handle)$SNP) expect_identical(.panelVariantFilter(handle, ids), ids) }) test_that(".panelVariantFilter drops panel-rare variants", { data(qtlDatasetExample) - handle <- qtlDatasetExample@genotypes + handle <- getGenotypeHandle(qtlDatasetExample) ids <- normalizeVariantId(getSnpInfo(handle)$SNP) loose <- .panelVariantFilter(handle, ids, mafCutoff = 0.05) tight <- .panelVariantFilter(handle, ids, mafCutoff = 0.2) @@ -4483,7 +4483,7 @@ test_that(".panelVariantFilter drops panel-rare variants", { test_that(".panelVariantFilter treats MAC as a MAF equivalent", { data(qtlDatasetExample) - handle <- qtlDatasetExample@genotypes + handle <- getGenotypeHandle(qtlDatasetExample) ids <- normalizeVariantId(getSnpInfo(handle)$SNP) nSamp <- getNSamples(handle) # macCutoff / (2 * nSamples) is the same threshold as mafCutoff. @@ -4499,7 +4499,7 @@ test_that(".panelVariantFilter treats MAC as a MAF equivalent", { test_that(".panelVariantFilter drops high-missingness variants", { data(qtlDatasetExample) - handle <- qtlDatasetExample@genotypes + handle <- getGenotypeHandle(qtlDatasetExample) ids <- normalizeVariantId(getSnpInfo(handle)$SNP) strict <- .panelVariantFilter(handle, ids, imissCutoff = 0) expect_lt(length(strict), length(ids)) @@ -4512,7 +4512,7 @@ test_that(".panelVariantFilter passes through ids absent from the panel", { # .ldFromSketch's `onMissing`; deciding it here too would let the two # disagree about the same variant. data(qtlDatasetExample) - handle <- qtlDatasetExample@genotypes + handle <- getGenotypeHandle(qtlDatasetExample) ids <- normalizeVariantId(getSnpInfo(handle)$SNP)[1:3] withGhost <- c("chr9:999:A:G", ids) expect_true(is_in( @@ -4523,7 +4523,7 @@ test_that(".panelVariantFilter passes through ids absent from the panel", { test_that(".panelVariantFilter handles empty and NULL input", { data(qtlDatasetExample) - handle <- qtlDatasetExample@genotypes + handle <- getGenotypeHandle(qtlDatasetExample) expect_length( .panelVariantFilter(handle, character(0), mafCutoff = 0.1), 0L diff --git a/tests/testthat/test_sumstatsQc.R b/tests/testthat/test_sumstatsQc.R index 85f5bea0..f06345e1 100644 --- a/tests/testthat/test_sumstatsQc.R +++ b/tests/testthat/test_sumstatsQc.R @@ -5672,204 +5672,6 @@ test_that("sliding_window_loop errors on infinite loop", { context("univariate_rss_diagnostics") -# A single-row QtlFineMappingResult standing in for the retired -# FineMappingRow. topLoci defaults to one row per variant: the row-payload -# builder requires the two to be aligned, where the entry tolerated an empty -# table beside a non-empty variant list. -.testFineMappingRow <- function( - variantIds, - susieFit = list(), - topLoci = NULL -) { - if (is.null(topLoci)) { - topLoci <- data.frame( - variant_id = variantIds, - pip = rep(0, length(variantIds)), - stringsAsFactors = FALSE - ) - } - QtlFineMappingResult( - study = "s1", - context = "c1", - trait = "t1", - method = "susie", - entry = list(fineMappingRow( - variantIds = variantIds, - susieFit = susieFit, - topLoci = topLoci - )) - ) -} - -# =========================================================================== -# getSusieResult -# =========================================================================== - -test_that("getSusieResult returns NULL for empty input", { - result <- getSusieResult(list()) - expect_null(result) -}) - -test_that("getSusieResult returns NULL when finemappingEntry missing", { - result <- getSusieResult(list(some_data = 42)) - expect_null(result) -}) - -test_that("getSusieResult returns trimmed result when present", { - mock_result <- list(pip = c(0.1, 0.5, 0.3), sets = list(cs = list())) - con_data <- list( - finemappingEntry = .testFineMappingRow( - variantIds = c("chr1:100:A:G", "chr1:200:C:T", "chr1:300:G:A"), - susieFit = mock_result - ) - ) - result <- getSusieResult(con_data) - expect_equal(result, mock_result) -}) - -# =========================================================================== -# extractTopPipInfo -# =========================================================================== - -test_that("extractTopPipInfo finds top PIP variant", { - con_data <- list( - finemappingEntry = .testFineMappingRow( - variantIds = c("chr1:100:A:G", "chr1:200:C:T", "chr1:300:G:A"), - susieFit = list(pip = c(0.1, 0.7, 0.2)) - ), - sumstats = list(z = c(1.0, 3.5, -0.5)) - ) - result <- extractTopPipInfo(con_data$finemappingEntry, con_data$sumstats) - expect_equal(result$top_variant, "chr1:200:C:T") - expect_equal(result$top_pip, 0.7) - expect_equal(result$top_z, 3.5) - expect_equal(result$top_variant_index, 2) - expect_true(is.na(result$cs_name)) - expect_true(is.na(result$variants_per_cs)) -}) - -test_that("extractTopPipInfo computes p_value from z", { - con_data <- list( - finemappingEntry = .testFineMappingRow( - variantIds = c("chr1:100:A:G", "chr1:200:C:T", "chr1:300:G:A"), - susieFit = list(pip = c(0.9, 0.05, 0.05)) - ), - sumstats = list(z = c(5.0, 0.5, -0.3)) - ) - result <- extractTopPipInfo(con_data$finemappingEntry, con_data$sumstats) - expected_pval <- pecotmr:::.zToPvalue(5.0) - expect_equal(result$p_value, expected_pval) -}) - -test_that("extractTopPipInfo handles ties by taking first max", { - con_data <- list( - finemappingEntry = .testFineMappingRow( - variantIds = c("chr1:100:A:G", "chr1:200:C:T", "chr1:300:G:A"), - susieFit = list(pip = c(0.5, 0.5, 0.5)) - ), - sumstats = list(z = c(1.0, 2.0, 3.0)) - ) - result <- extractTopPipInfo(con_data$finemappingEntry, con_data$sumstats) - expect_equal(result$top_variant_index, 1) - expect_equal(result$top_pip, 0.5) -}) - -# =========================================================================== -# extractCsInfo -# =========================================================================== - -test_that("extractCsInfo extracts single CS correctly", { - data(qtlSumStatsExample) - fe <- .testFineMappingRow( - variantIds = c("chr1:100:A:G", "chr1:200:C:T", "chr1:300:G:A"), - susieFit = list(sets = list(cs = list(L_1 = c(1, 2)))) - ) - top_loci_table <- data.frame( - variant_id = c("chr1:100:A:G", "chr1:200:C:T"), - pip = c(0.3, 0.8), - z = c(2.0, 4.5), - stringsAsFactors = FALSE - ) - # A single CS short-circuits (no between-CS correlation), so the unrelated - # ldSource is not consulted. - result <- extractCsInfo( - fe, - csNames = "L_1", - topLociTable = top_loci_table, - ldSource = qtlSumStatsExample - ) - expect_equal(nrow(result), 1) - expect_equal(result$cs_name, "L_1") - expect_equal(result$top_variant, "chr1:200:C:T") - expect_equal(result$top_pip, 0.8) - expect_equal(result$variants_per_cs, 2) - expect_true(is.na(result$cs_corr_max)) - expect_true(is.na(result$cs_corr_min)) - expect_false("cs_corr_1" %in% colnames(result)) -}) - -test_that("extractCsInfo builds correlation columns from the ldSource", { - data(qtlSumStatsExample) - ss <- qtlSumStatsExample - vids <- rownames(getLdSketch(ss)) - set.seed(1) - fit <- list( - sets = list(cs = list(L_1 = c(1L, 2L, 3L), L_2 = c(90L, 91L))), - pip = runif(length(vids)) - ) - tl <- data.frame( - variant_id = vids, - pip = fit$pip, - z = rnorm(length(vids)), - stringsAsFactors = FALSE - ) - fe <- .testFineMappingRow( - variantIds = vids, - susieFit = fit, - topLoci = tl - ) - result <- extractCsInfo( - fe, - csNames = c("L_1", "L_2"), - topLociTable = tl, - ldSource = ss - ) - expect_equal(nrow(result), 2) - expect_true(all( - c("cs_corr_1", "cs_corr_2", "cs_corr_max", "cs_corr_min") %in% - colnames(result) - )) - # The cs_corr_j columns are the columns of the computed between-CS matrix - # (symmetric; diagonal == 1), reduced on demand from the ldSource. - cc <- computeCsCorrelation(fe, ss) - expect_equal(result$cs_corr_1, unname(cc[, 1])) - expect_equal(result$cs_corr_2, unname(cc[, 2])) - expect_equal(result$cs_corr_max, rep(abs(cc[1, 2]), 2)) - expect_equal(result$cs_corr_min, rep(abs(cc[1, 2]), 2)) -}) - -test_that("extractCsInfo computes p_value from z-score", { - data(qtlSumStatsExample) - fe <- .testFineMappingRow( - variantIds = c("chr1:100:A:G", "chr1:200:C:T"), - susieFit = list(sets = list(cs = list(L_1 = c(1, 2)))) - ) - top_loci_table <- data.frame( - variant_id = c("chr1:100:A:G", "chr1:200:C:T"), - pip = c(0.9, 0.1), - z = c(5.0, 0.5), - stringsAsFactors = FALSE - ) - result <- extractCsInfo( - fe, - csNames = "L_1", - topLociTable = top_loci_table, - ldSource = qtlSumStatsExample - ) - expected_pval <- pecotmr:::.zToPvalue(5.0) - expect_equal(result$p_value, expected_pval, tolerance = 1e-10) -}) - # =========================================================================== # autoDecision # =========================================================================== @@ -6307,27 +6109,6 @@ test_that("slalom coerces a non-matrix X (data.frame) to a matrix", { expect_equal(nrow(result$data), n_snps) }) -# =========================================================================== -# getSusieResult: trimmed fit is empty -# =========================================================================== - -test_that("getSusieResult returns NULL when the trimmed susie fit is empty", { - conData <- list( - # An empty fit with variants present: topLoci must still be aligned - # row-for-row, which the entry did not enforce. - finemappingEntry = fineMappingRow( - variantIds = c("chr1:100:A:G", "chr1:200:C:T"), - susieFit = list(), - topLoci = data.frame( - variant_id = c("chr1:100:A:G", "chr1:200:C:T"), - pip = c(0, 0), - stringsAsFactors = FALSE - ) - ) - ) - expect_null(getSusieResult(conData)) -}) - # =========================================================================== # autoDecision: high-correlation tagging branch is reached # =========================================================================== diff --git a/tests/testthat/test_tupleSelectors.R b/tests/testthat/test_tupleSelectors.R index 4a46435a..2ab70715 100644 --- a/tests/testthat/test_tupleSelectors.R +++ b/tests/testthat/test_tupleSelectors.R @@ -1,20 +1,38 @@ context("tupleSelectors (internal row-selector helpers)") -# These helpers (.matchTupleRows, .tupleSelectRow, .tupleSelectRowGwasFmr) -# work on anything with `nrow(x)` and `x[[col]]` semantics, so the tests -# below use plain base-R data.frames rather than building S4 collections. +# These helpers read identity columns off `mcols(x)`, so the tests below use +# a real collection rather than a plain data.frame. `.ts_coll()` builds the +# smallest one that exercises them: one empty GRanges per row, identity +# columns in mcols, no payload class. + +setClass("TsTestCollection", contains = "RangedTupleList") + +.ts_coll <- function(..., stringsAsFactors = FALSE) { + cols <- list(...) + n <- if (length(cols) == 0L) 0L else length(cols[[1L]]) + grl <- GenomicRanges::GRangesList( + rep(list(GenomicRanges::GRanges()), n) + ) + if (length(cols) > 0L) { + S4Vectors::mcols(grl) <- do.call( + S4Vectors::DataFrame, + c(cols, list(check.names = FALSE)) + ) + } + methods::new("TsTestCollection", grl) +} # =========================================================================== # .matchTupleRows # =========================================================================== test_that(".matchTupleRows: empty keys returns every row index", { - df <- data.frame(study = c("s1", "s2"), method = c("susie", "lasso")) + df <- .ts_coll(study = c("s1", "s2"), method = c("susie", "lasso")) expect_equal(pecotmr:::.matchTupleRows(df, list()), c(1L, 2L)) }) test_that(".matchTupleRows: AND-matches across multiple (column, value) pairs", { - df <- data.frame( + df <- .ts_coll( study = c("s1", "s1", "s2"), context = c("c1", "c2", "c1"), stringsAsFactors = FALSE @@ -35,7 +53,7 @@ test_that(".matchTupleRows: AND-matches across multiple (column, value) pairs", # =========================================================================== test_that(".tupleSelectRow: zero-row input errors with the class label", { - empty <- data.frame( + empty <- .ts_coll( study = character(0), context = character(0), trait = character(0), @@ -56,7 +74,7 @@ test_that(".tupleSelectRow: zero-row input errors with the class label", { }) test_that(".tupleSelectRow: single-row collection returns 1L without selectors", { - one <- data.frame( + one <- .ts_coll( study = "s1", context = "c1", trait = "t1", @@ -67,7 +85,7 @@ test_that(".tupleSelectRow: single-row collection returns 1L without selectors", }) test_that(".tupleSelectRow: multi-row + missing selectors errors with row count", { - multi <- data.frame( + multi <- .ts_coll( study = c("s1", "s1"), context = c("c1", "c2"), trait = c("t1", "t1"), @@ -81,7 +99,7 @@ test_that(".tupleSelectRow: multi-row + missing selectors errors with row count" }) test_that(".tupleSelectRow: non-scalar selectors error", { - multi <- data.frame( + multi <- .ts_coll( study = c("s1", "s2"), context = c("c1", "c2"), trait = c("t1", "t2"), @@ -101,7 +119,7 @@ test_that(".tupleSelectRow: non-scalar selectors error", { }) test_that(".tupleSelectRow: matching tuple returns first row index", { - multi <- data.frame( + multi <- .ts_coll( study = c("s1", "s1"), context = c("c1", "c2"), trait = c("t1", "t1"), @@ -121,7 +139,7 @@ test_that(".tupleSelectRow: matching tuple returns first row index", { }) test_that(".tupleSelectRow: missing tuple errors with the 4-tuple in the message", { - multi <- data.frame( + multi <- .ts_coll( study = c("s1", "s1"), context = c("c1", "c2"), trait = c("t1", "t1"), @@ -145,7 +163,7 @@ test_that(".tupleSelectRow: missing tuple errors with the 4-tuple in the message # =========================================================================== test_that(".tupleSelectRowGwasFmr: zero-row input errors", { - empty <- data.frame( + empty <- .ts_coll( study = character(0), method = character(0), blockId = character(0), @@ -158,7 +176,7 @@ test_that(".tupleSelectRowGwasFmr: zero-row input errors", { }) test_that(".tupleSelectRowGwasFmr: single-row collection returns 1L", { - one <- data.frame( + one <- .ts_coll( study = "g1", method = "susie", blockId = "region_1", @@ -168,7 +186,7 @@ test_that(".tupleSelectRowGwasFmr: single-row collection returns 1L", { }) test_that(".tupleSelectRowGwasFmr: missing selectors on multi-row errors", { - multi <- data.frame( + multi <- .ts_coll( study = c("g1", "g2"), method = c("susie", "susie"), blockId = c("region_1", "region_1"), @@ -181,7 +199,7 @@ test_that(".tupleSelectRowGwasFmr: missing selectors on multi-row errors", { }) test_that(".tupleSelectRowGwasFmr: non-scalar region errors", { - multi <- data.frame( + multi <- .ts_coll( study = c("g1", "g1"), method = c("susie", "susie"), blockId = c("r1", "r2"), @@ -200,7 +218,7 @@ test_that(".tupleSelectRowGwasFmr: non-scalar region errors", { test_that(".tupleSelectRowGwasFmr: region disambiguates per-block rows", { # Same (study, method) across two regions; region picks the right row. - multi <- data.frame( + multi <- .ts_coll( study = c("g1", "g1"), method = c("susie", "susie"), blockId = c("chr22_1_100", "chr22_500_600"), @@ -218,7 +236,7 @@ test_that(".tupleSelectRowGwasFmr: region disambiguates per-block rows", { }) test_that(".tupleSelectRowGwasFmr: missing tuple errors and includes region in message", { - one <- data.frame( + one <- .ts_coll( study = "g1", method = "susie", blockId = "r1", @@ -238,7 +256,7 @@ test_that(".tupleSelectRowGwasFmr: missing tuple errors and includes region in m test_that(".tupleSelectRowGwasFmr: ambiguous multi-match (no region) lists candidates", { # Two rows share (study, method); .tupleSelectRowGwasFmr should error # listing the available blockIds since the caller didn't disambiguate. - multi <- data.frame( + multi <- .ts_coll( study = c("g1", "g1"), method = c("susie", "susie"), blockId = c("region_A", "region_B"), @@ -255,7 +273,7 @@ test_that(".tupleSelectRowGwasFmr: ambiguous multi-match (no region) lists candi # =========================================================================== test_that(".fmrRowsMatching: no selectors returns every row", { - df <- data.frame( + df <- .ts_coll( study = c("s1", "s1"), context = c("c1", "c2"), trait = c("t1", "t1"), @@ -266,7 +284,7 @@ test_that(".fmrRowsMatching: no selectors returns every row", { }) test_that(".fmrRowsMatching: matches a subset without erroring on ambiguity", { - df <- data.frame( + df <- .ts_coll( study = c("s1", "s1", "s2"), context = c("c1", "c2", "c1"), trait = c("t1", "t1", "t1"), @@ -280,7 +298,7 @@ test_that(".fmrRowsMatching: matches a subset without erroring on ambiguity", { test_that(".fmrRowsMatching: selectors on absent columns are ignored", { # GWAS-shaped frame has no context/trait column; passing context must not # error and must not constrain the result. - gwas <- data.frame( + gwas <- .ts_coll( study = c("g1", "g1"), method = c("susie", "susie"), blockId = c("r1", "r2"), @@ -290,7 +308,7 @@ test_that(".fmrRowsMatching: selectors on absent columns are ignored", { }) test_that(".fmrRowsMatching: `region` matches the blockId column", { - gwas <- data.frame( + gwas <- .ts_coll( study = c("g1", "g1"), method = c("susie", "susie"), blockId = c("r1", "r2"), @@ -300,7 +318,7 @@ test_that(".fmrRowsMatching: `region` matches the blockId column", { }) test_that(".fmrRowsMatching: a vector selector matches any listed value", { - df <- data.frame( + df <- .ts_coll( study = c("s1", "s2", "s3"), context = c("c1", "c2", "c3"), trait = c("t1", "t1", "t1"), @@ -318,7 +336,7 @@ test_that(".fmrRowsMatching: a vector selector matches any listed value", { # =========================================================================== test_that(".fmrRowMetadata: emits all five identity columns, NA-filling absent ones", { - qtl <- data.frame( + qtl <- .ts_coll( study = c("s1", "s1"), context = c("c1", "c2"), trait = c("t1", "t1"), @@ -335,7 +353,7 @@ test_that(".fmrRowMetadata: emits all five identity columns, NA-filling absent o }) test_that(".fmrRowMetadata: GWAS frame NA-fills context/trait, keeps blockId", { - gwas <- data.frame( + gwas <- .ts_coll( study = c("g1", "g1"), method = c("susie", "susie"), blockId = c("r1", "r2"), @@ -348,7 +366,7 @@ test_that(".fmrRowMetadata: GWAS frame NA-fills context/trait, keeps blockId", { }) test_that(".fmrRowMetadata: zero-row input yields a zero-row 5-column frame", { - empty <- data.frame( + empty <- .ts_coll( study = character(0), method = character(0), blockId = character(0), @@ -401,8 +419,8 @@ test_that(".rbindCollections: all-NULL input returns NULL", { expect_null(pecotmr:::.rbindCollections(list(NULL, NULL))) }) -test_that(".getRegionColumn: an absent region column yields an empty GRanges", { - gr <- pecotmr:::.getRegionColumn(S4Vectors::DataFrame(a = 1:2)) +test_that(".getRegionColumn: a zero-row collection yields an empty GRanges", { + gr <- pecotmr:::.getRegionColumn(.ts_coll(study = character(0))) expect_s4_class(gr, "GRanges") expect_length(gr, 0L) }) @@ -453,13 +471,14 @@ test_that(".appendTraitPosCol: traitPos must be a GRanges of matching length", { test_that(".validateTraitPosColumn: reports non-GRanges and wrong-length traitPos", { expect_equal( - pecotmr:::.validateTraitPosColumn(S4Vectors::DataFrame( - traitPos = c("x", "y") - )), + pecotmr:::.validateTraitPosColumn(.ts_coll(traitPos = c("x", "y"))), "'traitPos' column must be a GRanges" ) - bad <- S4Vectors::DataFrame(a = 1:2) - bad@listData$traitPos <- GenomicRanges::GRanges( + # A one-range traitPos beside two rows: assigned through the mcols + # listData because the parallel-length check would reject it otherwise, + # which is exactly the state the validator has to catch. + bad <- .ts_coll(a = 1:2) + S4Vectors::mcols(bad)@listData$traitPos <- GenomicRanges::GRanges( "chr1", IRanges::IRanges(1, 1) ) @@ -468,7 +487,7 @@ test_that(".validateTraitPosColumn: reports non-GRanges and wrong-length traitPo "'traitPos' column must have one range per row" ) expect_length( - pecotmr:::.validateTraitPosColumn(S4Vectors::DataFrame(a = 1L)), + pecotmr:::.validateTraitPosColumn(.ts_coll(a = 1L)), 0L ) })