diff --git a/NAMESPACE b/NAMESPACE index 4d7e86aa..90d26e3f 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -134,8 +134,11 @@ export(getGenome) export(getGenotypeCovariates) export(getGenotypes) export(getH2) +export(getH2Se) export(getImissCutoff) export(getInSample) +export(getIntercept) +export(getInterceptSe) export(getKeepIndel) export(getKeepVariants) export(getLbf) @@ -153,6 +156,7 @@ export(getMethodNames) export(getMixtureWeights) export(getN) export(getNRef) +export(getNSnps) export(getP) export(getPhenotypeCovariates) export(getPhenotypes) @@ -181,6 +185,7 @@ export(getSusieFit) export(getSusieResult) export(getTauBlocks) export(getTopLoci) +export(getTraitName) export(getTraitNames) export(getTraitPosition) export(getTraitRun) @@ -365,8 +370,11 @@ exportMethods(getGenome) exportMethods(getGenotypeCovariates) exportMethods(getGenotypes) exportMethods(getH2) +exportMethods(getH2Se) exportMethods(getImissCutoff) exportMethods(getInSample) +exportMethods(getIntercept) +exportMethods(getInterceptSe) exportMethods(getKeepIndel) exportMethods(getKeepVariants) exportMethods(getLbf) @@ -384,6 +392,7 @@ exportMethods(getMethodNames) exportMethods(getMixtureWeights) exportMethods(getN) exportMethods(getNRef) +exportMethods(getNSnps) exportMethods(getP) exportMethods(getPhenotypeCovariates) exportMethods(getPhenotypes) @@ -402,7 +411,6 @@ exportMethods(getScoreStats) exportMethods(getSe) exportMethods(getSignificantQtls) exportMethods(getSnpIdx) -exportMethods(getSnpInfo) exportMethods(getStandardized) exportMethods(getStudy) exportMethods(getSumStats) @@ -411,6 +419,7 @@ exportMethods(getSusieAlpha) exportMethods(getSusieFit) exportMethods(getTauBlocks) exportMethods(getTopLoci) +exportMethods(getTraitName) exportMethods(getTraitNames) exportMethods(getTraitPosition) exportMethods(getTraitRun) @@ -552,6 +561,7 @@ importFrom(purrr, map_dfr, map_int, map_lgl, + partial, reduce, set_names ) diff --git a/R/AllGenerics.R b/R/AllGenerics.R index 2742e310..45e9db69 100644 --- a/R/AllGenerics.R +++ b/R/AllGenerics.R @@ -1162,9 +1162,14 @@ setGeneric("getFullFit", function(x, ...) standardGeneric("getFullFit")) setGeneric("getCvFits", function(x, ...) standardGeneric("getCvFits")) #' @title Get Method Names -#' @description Extract method names from a collection class. -#' @param x A \code{FineMappingResult} or \code{TwasWeights} object. +#' @description Extract the method name(s) an object was produced with: one +#' per entry for a collection, a single name for one estimate. +#' @param x A \code{FineMappingResult}, \code{TwasWeights}, +#' \code{CtwasResult} or \code{H2Estimate} object. #' @return Character vector. +#' @examples +#' data(h2EstimateExample) +#' getMethodNames(h2EstimateExample) #' @export setGeneric("getMethodNames", function(x) standardGeneric("getMethodNames")) @@ -1441,9 +1446,16 @@ setGeneric("getScaleResiduals", function(x) { #' @title Get SNP Info #' @description Return the cached SNP metadata data.frame (columns: SNP, CHR, -#' BP, A1, A2, optionally MAF). +#' BP, A1, A2, fileIdx, optionally MAF). #' @param x An object carrying cached SNP metadata. #' @return A data.frame. +#' @details The genotype handle is the seed layer behind a panel and +#' is not part of the public interface; obtain a panel from +#' \code{readGenotypes()} and ask it directly. This table is the READ-PATH +#' view: unlike the panel's \code{rowRanges()} it carries \code{fileIdx}, +#' the on-disk variant position that keeps reads correct after a handle has +#' been row-subset. Callers that only need chrom/pos/alleles should use +#' \code{rowRanges()} on the panel instead. #' @keywords internal setGeneric("getSnpInfo", function(x) standardGeneric("getSnpInfo")) @@ -1708,6 +1720,61 @@ setGeneric("getTauBlocks", function(x) standardGeneric("getTauBlocks")) #' @export setGeneric("getH2", function(x) standardGeneric("getH2")) +#' @title Get Heritability Standard Error +#' @description Return the standard error of the global SNP heritability +#' estimate carried by an \code{H2Estimate}. +#' @param x An \code{H2Estimate}. +#' @return Numeric (length 1). +#' @examples +#' data(h2EstimateExample) +#' getH2Se(h2EstimateExample) +#' @export +setGeneric("getH2Se", function(x) standardGeneric("getH2Se")) + +#' @title Get LD-Score Regression Intercept +#' @description Return the regression intercept carried by an +#' \code{H2Estimate}. Values above 1 indicate confounding or sample overlap +#' rather than polygenic signal. +#' @param x An \code{H2Estimate}. +#' @return Numeric (length 1). +#' @examples +#' data(h2EstimateExample) +#' getIntercept(h2EstimateExample) +#' @export +setGeneric("getIntercept", function(x) standardGeneric("getIntercept")) + +#' @title Get Intercept Standard Error +#' @description Return the standard error of the regression intercept carried +#' by an \code{H2Estimate}. +#' @param x An \code{H2Estimate}. +#' @return Numeric (length 1). +#' @examples +#' data(h2EstimateExample) +#' getInterceptSe(h2EstimateExample) +#' @export +setGeneric("getInterceptSe", function(x) standardGeneric("getInterceptSe")) + +#' @title Get Variant Count +#' @description Return the number of variants the heritability estimate was +#' computed over. +#' @param x An \code{H2Estimate}. +#' @return Integer (length 1). +#' @examples +#' data(h2EstimateExample) +#' getNSnps(h2EstimateExample) +#' @export +setGeneric("getNSnps", function(x) standardGeneric("getNSnps")) + +#' @title Get Trait Name +#' @description Return the trait label carried by an \code{H2Estimate}. +#' @param x An \code{H2Estimate}. +#' @return Character (length 1). +#' @examples +#' data(h2EstimateExample) +#' getTraitName(h2EstimateExample) +#' @export +setGeneric("getTraitName", function(x) standardGeneric("getTraitName")) + # Internal generics for the unified joint-analysis engine (see R/JointGroup.R # and dev/jointSpecification-s4-refactor.md). Not exported: the engine and its # fitters are package-internal machinery. diff --git a/R/GenotypeHandle.R b/R/GenotypeHandle.R index 5d6369bc..668604b7 100644 --- a/R/GenotypeHandle.R +++ b/R/GenotypeHandle.R @@ -680,7 +680,7 @@ GenotypeHandle <- function( } #' @rdname getSnpInfo -#' @export +#' @keywords internal setMethod("getSnpInfo", "GenotypeHandle", function(x) x@snpInfo) #' @rdname getFormat diff --git a/R/H2Estimate.R b/R/H2Estimate.R index 378b461f..8ea3e5ee 100644 --- a/R/H2Estimate.R +++ b/R/H2Estimate.R @@ -64,6 +64,30 @@ setClass( #' @export setMethod("getH2", "H2Estimate", function(x) x@h2) +#' @rdname getH2Se +#' @export +setMethod("getH2Se", "H2Estimate", function(x) x@h2Se) + +#' @rdname getIntercept +#' @export +setMethod("getIntercept", "H2Estimate", function(x) x@intercept) + +#' @rdname getInterceptSe +#' @export +setMethod("getInterceptSe", "H2Estimate", function(x) x@interceptSe) + +#' @rdname getNSnps +#' @export +setMethod("getNSnps", "H2Estimate", function(x) x@nSnps) + +#' @rdname getTraitName +#' @export +setMethod("getTraitName", "H2Estimate", function(x) x@traitName) + +#' @rdname getMethodNames +#' @export +setMethod("getMethodNames", "H2Estimate", function(x) x@method) + #' @rdname getTauBlocks #' @export setMethod("getTauBlocks", "H2Estimate", function(x) x@tauBlocks) diff --git a/R/JointGroup.R b/R/JointGroup.R index 87a116cc..4c2ef6de 100644 --- a/R/JointGroup.R +++ b/R/JointGroup.R @@ -65,7 +65,7 @@ setClass( setClass( "IndividualJointGroup", contains = "JointGroup", - representation(X = "matrix", Y = "matrix", pos = "numeric"), + representation(X = "matrix", Y = "matrix", traitPos = "numeric"), validity = function(object) { errors <- character() if (nrow(object@X) != nrow(object@Y)) { @@ -74,10 +74,16 @@ setClass( if (ncol(object@Y) != nrow(object@conditions)) { errors <- c(errors, "ncol(Y) must equal nrow(conditions)") } - if (length(object@pos) > 0L && length(object@pos) != ncol(object@Y)) { + if ( + length(object@traitPos) > 0L && + length(object@traitPos) != ncol(object@Y) + ) { errors <- c( errors, - "when set, 'pos' must have one entry per Y column" + str_c( + "when set, 'traitPos' must have one entry per TRAIT ", + "(Y column), not per variant" + ) ) } if (length(errors) == 0L) TRUE else errors @@ -88,21 +94,26 @@ setClass( setClass( "SumStatsJointGroup", contains = "JointGroup", - representation(Z = "matrix", R = "matrix", N = "numeric"), + # The LD REFERENCE, not the LD matrix: R is variants x variants, and + # every group in an enumeration would hold its own dense copy long before + # any of them is fitted. The sketch is a lazy panel, so the matrix is + # derived once per group at fit time -- which is what the univariate RSS + # path in fineMappingPipeline already does. + representation(Z = "matrix", ldSketch = "LdSketchOrNULL", N = "numeric"), validity = function(object) { errors <- character() - if (nrow(object@R) != ncol(object@R)) { - errors <- c(errors, "'R' (LD) must be square") + if (ncol(object@Z) != nrow(object@conditions)) { + errors <- c(errors, "ncol(Z) must equal nrow(conditions)") } - if (nrow(object@Z) != nrow(object@R)) { + if (is.null(rownames(object@Z))) { errors <- c( errors, - "'Z' rows (variants) must match the 'R' dimension" + str_c( + "'Z' must carry the variant ids as rownames: they are ", + "what the LD matrix is derived over" + ) ) } - if (ncol(object@Z) != nrow(object@conditions)) { - errors <- c(errors, "ncol(Z) must equal nrow(conditions)") - } if (length(errors) == 0L) TRUE else errors } ) @@ -159,14 +170,36 @@ setClass("TwasJointPipeline", contains = "JointPipeline") # @noRd .jgY <- function(g) g@Y +# The genomic midpoint of each TRAIT, one per Y column -- fsusie's functional +# domain coordinate, not a variant position. Named apart from `pos`, which +# everywhere else in the package means the position of a variant. # @noRd -.jgPos <- function(g) g@pos +.jgTraitPos <- function(g) g@traitPos # @noRd .jgZ <- function(g) g@Z # @noRd -.jgR <- function(g) g@R +.jgLdSketch <- function(g) g@ldSketch + +# The variants a summary-statistics group covers, in Z's row order -- which is +# the order the derived LD matrix comes back in. +# @noRd +.jgVariantIds <- function(g) rownames(.jgZ(g)) + +# The group's LD matrix, derived from its sketch. +# +# Callers must derive ONCE per group and pass the result down: the per- +# condition entry builder runs once per Z column, so calling this from there +# would recompute the same matrix for every condition. +# @noRd +.jgLdMatrix <- function(g) { + .ldFromSketch( + .jgLdSketch(g), + .jgVariantIds(g), + label = "jointEngine" + ) +} # @noRd .jgN <- function(g) g@N diff --git a/R/LdData.R b/R/LdData.R index 6eae9ac3..bcefd808 100644 --- a/R/LdData.R +++ b/R/LdData.R @@ -15,6 +15,14 @@ NULL #' correlation matrix or a \code{GenotypeHandle} (or list of handles for #' mixture panels) for lazy genotype/correlation access. #' +#' An \code{LdData} \strong{is} a \code{GRanges} over the variants it +#' covers -- carrying A1, A2, variant_id and optionally allele_freq, +#' variance and n_nomiss as metadata columns -- so \code{length()}, +#' \code{seqnames()} and \code{start()} answer directly, as they do for +#' \code{\link{LdScore}} and \code{\link{LdEigen}}. The correlation may +#' be NULL, in which case those ranges are the object's only record of which +#' variants it describes. +#' #' @slot correlation A correlation matrix, a list of per-block matrices #' (block-diagonal LD), or NULL if genotypes are available and R should be #' computed on demand. @@ -26,8 +34,6 @@ NULL #' @slot snpIdx Integer vector of 1-based SNP indices into the handle's #' \code{snpInfo}. NULL when correlation is pre-computed, or when the #' source is a matrix (which is already the subset). -#' @slot variants A \code{GRanges} object with variant metadata (A1, A2, -#' variant_id, and optionally allele_freq, variance, n_nomiss). #' @slot blockMetadata A \code{GRanges} of blocks or a \code{data.frame} with #' block boundary information. #' @slot nRef Integer, reference panel sample size. @@ -40,11 +46,11 @@ NULL #' @export setClass( "LdData", + contains = "GRanges", representation( correlation = "LdCorrelation", genotypeHandle = "LdGenotypeSource", snpIdx = "LdSnpIndex", - variants = "GRanges", blockMetadata = "LdBlockMetadata", nRef = "integer", mixtureWeights = "LdMixtureWeights" @@ -60,8 +66,8 @@ setClass( ) ) } - if (length(object@variants) == 0) { - errors <- c(errors, "'variants' must not be empty") + if (length(object) == 0) { + errors <- c(errors, "an LdData must cover >= 1 variant") } errors <- c(errors, .ldCheckGenotypeSource(object@genotypeHandle)) errors <- c(errors, .ldCheckCorrelation(object@correlation)) @@ -98,10 +104,29 @@ setClass( } ) +#' @describeIn LdData-class Refused. Subsetting would narrow the variants +#' while \code{correlation} -- variant-by-variant, and \code{snpIdx}, which +#' indexes the reference panel -- stayed as they were, leaving an LD matrix +#' describing variants the object no longer has. Build the LdData over the +#' variant set you want instead. +#' @param x An \code{LdData}. +#' @param i,j,... Subscripts; any use is an error. +#' @param drop Ignored. +#' @return Nothing: this method always signals an error. +#' @export +setMethod("[", "LdData", function(x, i, j, ..., drop = TRUE) { + abort(glue( + "an LdData cannot be subset: its correlation is variant-by-variant ", + "and its snpIdx addresses the reference panel, so narrowing the ", + "variants would leave both describing variants that are no longer ", + "present. Build the LdData over the variant set you want instead." + )) +}) + #' @rdname show-methods #' @export setMethod("show", "LdData", function(object) { - n_var <- length(object@variants) + n_var <- length(object) has_R <- !is.null(object@correlation) has_geno <- !is.null(object@genotypeHandle) r_type <- if (has_R && is.list(object@correlation)) { @@ -201,10 +226,10 @@ LdData <- function( ) { obj <- new( "LdData", + variants, correlation = correlation, genotypeHandle = .ldDataGenotypeSource(genotypeHandle), snpIdx = if (is.null(snpIdx)) NULL else as.integer(snpIdx), - variants = variants, blockMetadata = blockMetadata, nRef = as.integer(nRef), mixtureWeights = mixtureWeights @@ -214,8 +239,7 @@ LdData <- function( } # Internal: convert a refPanel data.frame (chrom/pos/A1/A2/variant_id, with -# optional allele_freq/variance/n_nomiss) into the GRanges form used by the -# LdData `variants` slot. +# optional allele_freq/variance/n_nomiss) into the GRanges an LdData is. .refPanelToGranges <- function(refPanel) { chr <- withChrPrefix(refPanel$chrom) pos <- as.integer(refPanel$pos) @@ -303,13 +327,15 @@ setMethod("hasGenotypes", "LdData", function(x) { #' @rdname getVariantIds #' @export setMethod("getVariantIds", "LdData", function(x, ...) { - mcols(x@variants)$variant_id + mcols(x)$variant_id }) #' @rdname getVariantInfo #' @export setMethod("getVariantInfo", "LdData", function(x) { - x@variants + # A plain GRanges view: callers subset and re-wrap it, and carrying the + # LD payload along would make those copies quietly expensive. + as(x, "GRanges") }) #' @rdname getBlockMetadata @@ -321,9 +347,9 @@ setMethod("getBlockMetadata", "LdData", function(x) { #' @rdname getRefPanel #' @export setMethod("getRefPanel", "LdData", function(x) { - mc <- as_tibble(as.data.frame(mcols(x@variants))) - mc$chrom <- as.character(seqnames(x@variants)) - mc$pos <- start(x@variants) + mc <- as_tibble(as.data.frame(mcols(x))) + mc$chrom <- as.character(seqnames(x)) + mc$pos <- start(x) mc }) diff --git a/R/causalInferencePipeline.R b/R/causalInferencePipeline.R index cdde3596..1973e60a 100644 --- a/R/causalInferencePipeline.R +++ b/R/causalInferencePipeline.R @@ -733,7 +733,7 @@ causalInferencePipeline <- function( wSub <- weights[m$idxA] * m$sign zSub <- gwasDf$z[m$idxB] gwasIds <- gwasDf$variant_id[m$idxB] - ldMat <- .cipLdFromSketch(gwasLd, gwasIds) + ldMat <- .ldFromSketch(gwasLd, gwasIds, label = "causalInferencePipeline") res <- twasZ(weights = wSub, z = zSub, R = ldMat) zMat <- res$Z @@ -777,13 +777,6 @@ causalInferencePipeline <- function( warn(msg) } -# Build an LD correlation matrix for a given variant subset from an -# LD sketch. Thin wrapper over the shared `.ldFromSketch` helper -# (R/ld.R). -.cipLdFromSketch <- function(ldSketch, variantIds) { - .ldFromSketch(ldSketch, variantIds, label = "causalInferencePipeline") -} - # Compute the Wald-ratio IVW MR estimate for a single tuple. Uses the # FineMappingRow's topLoci as the instrumental variable source: each # variant with PIP > pipCutoff contributes one ratio = beta_y / beta_x. diff --git a/R/crossValidation.R b/R/crossValidation.R index 2a589485..2e85d295 100644 --- a/R/crossValidation.R +++ b/R/crossValidation.R @@ -441,7 +441,7 @@ if (verbose < 1) { return(invisible(NULL)) } - s <- suppressWarnings(stats::sd(pred, na.rm = TRUE)) + s <- stats::sd(pred, na.rm = TRUE) if (is.na(s) || s == 0) { msg <- glue( "Predicted values for condition {r} using {mk} have zero ", diff --git a/R/ctwasPipeline.R b/R/ctwasPipeline.R index 49140632..9e59d934 100644 --- a/R/ctwasPipeline.R +++ b/R/ctwasPipeline.R @@ -2168,19 +2168,17 @@ asCtwasResult <- function(finemapResult, keepSnps = FALSE) { # supported by the downstream cTWAS model. # @noRd .ctwasSnpInfoForBlock <- function(gwasLd) { - snpInfo <- .ldSketchSnpInfo(gwasLd) - chr <- as.integer( - str_remove(as.character(snpInfo$CHR), regex("^chr", ignore_case = TRUE)) - ) + gr <- .ldSketchRanges(gwasLd) + mc <- S4Vectors::mcols(gr) # Base data.frame (not tibble): ctwas indexes snp_map positionally # (df[, "pos"] -> vector for region-bound min/max); a tibble column is a # 1-col list and errors min(). data.frame( - chrom = chr, - id = as.character(snpInfo$SNP), - pos = as.integer(snpInfo$BP), - alt = as.character(snpInfo$A1), - ref = as.character(snpInfo$A2), + chrom = as.integer(.ldSketchChrom(gwasLd)), + id = as.character(mc$SNP), + pos = as.integer(GenomicRanges::start(gr)), + alt = as.character(mc$A1), + ref = as.character(mc$A2), stringsAsFactors = FALSE ) } @@ -2198,8 +2196,8 @@ asCtwasResult <- function(finemapResult, keepSnps = FALSE) { # @noRd .ctwasComputeFullPanelLd <- function(gwasLd) { snpInfoCtwas <- .ctwasSnpInfoForBlock(gwasLd) - geno <- .dosageMatrix( - .ldSketchHandle(gwasLd), + geno <- .ldSketchDosage( + gwasLd, seq_len(nrow(snpInfoCtwas)), meanImpute = TRUE ) diff --git a/R/fineMappingPipeline.R b/R/fineMappingPipeline.R index e8e3fff1..aa498b3f 100644 --- a/R/fineMappingPipeline.R +++ b/R/fineMappingPipeline.R @@ -853,7 +853,11 @@ setGeneric("fineMappingPipeline", function(data, ...) { .fmGwasFitRows <- function(p, gr, zn, st, blockId, toRun) { z <- zn$z names(z) <- zn$variantIds - ldMat <- .fmLdFromSketch(p$ldSketch, zn$variantIds) + ldMat <- .ldFromSketch( + p$ldSketch, + zn$variantIds, + label = "fineMappingPipeline" + ) ents <- .fmFitRssBlock( z, ldMat, @@ -938,7 +942,11 @@ setGeneric("fineMappingPipeline", function(data, ...) { ) z <- zn$z names(z) <- zn$variantIds - ldMat <- .fmLdFromSketch(p$ldSketch, zn$variantIds) + ldMat <- .ldFromSketch( + p$ldSketch, + zn$variantIds, + label = "fineMappingPipeline" + ) ents <- .fmFitRssBlock( z, ldMat, @@ -1063,15 +1071,6 @@ combineFineMappingResults <- function(..., ldSketch = NULL) { } -# Build an LD correlation matrix from an LD sketch genotype handle for a -# specific variant subset. Thin wrapper over the shared `.ldFromSketch` -# helper (R/ld.R). -# @noRd -.fmLdFromSketch <- function(ldSketch, variantIds) { - .ldFromSketch(ldSketch, variantIds, label = ".fmLdFromSketch") -} - - # Wrap one finemapping fit into a FineMappingRow via the surviving # post-processing helpers (postprocessFinemappingFits + # formatFinemappingOutput). Returns a bare FineMappingRow payload, ready @@ -3186,7 +3185,7 @@ setMethod( is.null(rFinite) && (isTRUE(serFallback) || !identical(rMismatch, "none")) ) { - getNSamples(.ldSketchHandle(ldSketch)) + .ldSketchNSamples(ldSketch) } else { rFinite } diff --git a/R/fineMappingWrappers.R b/R/fineMappingWrappers.R index 4f90a638..a0ee4506 100644 --- a/R/fineMappingWrappers.R +++ b/R/fineMappingWrappers.R @@ -1425,7 +1425,8 @@ buildTopLoci <- function( # dropped its N to NA. For the list case fitN stays NA and the per-variant # `n` (threaded from the RSS effective sample size) fills the N column. hasOutcomeMatrix <- !is.null(dataY) && - (is.matrix(dataY) || is.data.frame(dataY) || + (is.matrix(dataY) || + is.data.frame(dataY) || (is.atomic(dataY) && is.numeric(dataY))) dataYMat <- if (hasOutcomeMatrix) as.matrix(dataY) else NULL fitN <- if (is.null(dataYMat)) NA_integer_ else as.integer(nrow(dataYMat)) diff --git a/R/genotypeIo.R b/R/genotypeIo.R index e4b1837e..f344c6cb 100644 --- a/R/genotypeIo.R +++ b/R/genotypeIo.R @@ -437,7 +437,8 @@ extractBlockGenotypes <- function(handle, snpIdx, meanImpute = TRUE) { as.matrix(genotypeCovariates) } SummarizedExperiment::colData(panel) <- .genotypeColData( - gCov, colnames(panel) + gCov, + colnames(panel) ) panel } @@ -729,10 +730,13 @@ extractBlockGenotypes <- function(handle, snpIdx, meanImpute = TRUE) { # layout. The shared form of the t(assay(extractBlockGenotypes(...))) idiom. # @noRd .dosageMatrix <- function(handle, snpIdx, meanImpute = TRUE) { - t(SummarizedExperiment::assay( - extractBlockGenotypes(handle, snpIdx, meanImpute = meanImpute), - "dosage" - )) + block <- extractBlockGenotypes(handle, snpIdx, meanImpute = meanImpute) + # NULL, not an error, when the reader declines the block: callers read + # that as "no genotypes here" (computeLd answers with an identity). + if (is.null(block)) { + return(NULL) + } + t(SummarizedExperiment::assay(block, "dosage")) } # Open a GDS read-only, guarantee it is closed on exit, and return fn(gds). diff --git a/R/jointEngine.R b/R/jointEngine.R index e6930fdb..3c14a321 100644 --- a/R/jointEngine.R +++ b/R/jointEngine.R @@ -256,7 +256,7 @@ setMethod( # cross-trait). One per-condition entry per trait, with an optional CV slice. # @noRd .jointFitFsusie <- function(group, Xc, Yc, nCond, cfg, args) { - if (length(.jgPos(group)) != nCond) { + if (length(.jgTraitPos(group)) != nCond) { msg <- glue( "fitJointGroup: fsusie requires per-trait positions ('pos'); ", "it is cross-trait individual-level only." @@ -265,7 +265,7 @@ setMethod( } verbose <- if (is.null(cfg$verbose)) 1 else cfg$verbose fitArgs <- .fmMergeUserArgs( - list(X = Xc, Y = Yc, pos = .jgPos(group)), + list(X = Xc, Y = Yc, pos = .jgTraitPos(group)), "fsusie", args$methodArgs[["fsusie"]] ) @@ -304,7 +304,7 @@ setMethod( cvFolds, samplePartition = cfg$samplePartition, coverage = cfg$coverage, - pos = .jgPos(group), + pos = .jgTraitPos(group), verbose = verbose, numThreads = if (is.null(cfg$cvThreads)) 1L else cfg$cvThreads, seed = cfg$seed @@ -535,13 +535,18 @@ setMethod( abort(msg) } cfg <- .jpConfig(pipeline) - fit <- .jointFitMvsusieRss(group, cfg, args) + # Derived ONCE here and threaded down: .jointRssEntry runs per + # condition, so deriving inside it would rebuild the same matrix for + # every column of Z. + ldMat <- .jgLdMatrix(group) + fit <- .jointFitMvsusieRss(group, ldMat, cfg, args) # One per-condition entry (RSS has no sample folds). map( seq_len(ncol(.jgZ(group))), .jointRssEntry, fit = fit, group = group, + ldMat = ldMat, cfg = cfg ) } @@ -550,7 +555,7 @@ setMethod( # mvSuSiE-RSS joint fit over summary statistics, with the data-driven reweighted # prior. Returns the class-tagged fit. # @noRd -.jointFitMvsusieRss <- function(group, cfg, args) { +.jointFitMvsusieRss <- function(group, ldMat, cfg, args) { ddCut <- if (is.null(cfg$dataDrivenPriorWeightsCutoff)) { 1e-10 } else { @@ -570,7 +575,7 @@ setMethod( ) mvBaseArgs <- list( Z = .jgZ(group), - R = .jgR(group), + R = ldMat, N = as.numeric(stats::median(.jgN(group))), prior_variance = mvPrior$priorVariance, coverage = cfg$coverage @@ -589,11 +594,11 @@ setMethod( # One RSS per-condition FineMappingRow (csInput = "Xcorr"). # @noRd -.jointRssEntry <- function(r, fit, group, cfg) { +.jointRssEntry <- function(r, fit, group, ldMat, cfg) { .fmPostprocessOne( fit = fit, method = "mvsusie", - dataX = .jgR(group), + dataX = ldMat, dataY = NULL, conditionIdx = r, coverage = cfg$coverage, @@ -958,7 +963,7 @@ setMethod( rfd <- if (is.null(cfg$retainFitDetail)) "slim" else cfg$retainFitDetail weights <- mrmashRssWeights( stat = list(z = .jgZ(group), n = .jgN(group)), - LD = .jgR(group), + LD = .jgLdMatrix(group), retainFit = TRUE, fitDetail = rfd ) @@ -1093,7 +1098,6 @@ setMethod("construct", "TwasJointPipeline", function(pipeline, records, ...) { ldSketch = ldSketch, cutoffs = args$cutoffs ) - ldMat <- .fmLdFromSketch(ldSketch, jz$variantIds) groups[[length(groups) + 1L]] <- new( "SumStatsJointGroup", conditions = tibble( @@ -1102,7 +1106,7 @@ setMethod("construct", "TwasJointPipeline", function(pipeline, records, ...) { trait = tid ), Z = jz$Z, - R = ldMat, + ldSketch = ldSketch, N = jz$nVec ) } @@ -1139,7 +1143,7 @@ setMethod("construct", "TwasJointPipeline", function(pipeline, records, ...) { # mvsusie ignores them. Matches the trait order of Y. rr <- SummarizedExperiment::rowRanges(xy$se) rr <- rr[match(colnames(xy$Y), rownames(xy$se))] - pos <- (GenomicRanges::start(rr) + GenomicRanges::end(rr)) / 2 + traitPos <- (GenomicRanges::start(rr) + GenomicRanges::end(rr)) / 2 groups[[length(groups) + 1L]] <- new( "IndividualJointGroup", conditions = tibble( @@ -1149,7 +1153,7 @@ setMethod("construct", "TwasJointPipeline", function(pipeline, records, ...) { ), X = xy$X, Y = xy$Y, - pos = as.numeric(pos) + traitPos = as.numeric(traitPos) ) } groups @@ -1181,7 +1185,6 @@ setMethod("construct", "TwasJointPipeline", function(pipeline, records, ...) { ldSketch = ldSketch, cutoffs = args$cutoffs ) - ldMat <- .fmLdFromSketch(ldSketch, jz$variantIds) groups[[length(groups) + 1L]] <- new( "SumStatsJointGroup", conditions = tibble( @@ -1190,7 +1193,7 @@ setMethod("construct", "TwasJointPipeline", function(pipeline, records, ...) { trait = trNames ), Z = jz$Z, - R = ldMat, + ldSketch = ldSketch, N = jz$nVec ) } @@ -1269,7 +1272,7 @@ setMethod("construct", "TwasJointPipeline", function(pipeline, records, ...) { trait = tid ), Z = jz$Z, - R = .fmLdFromSketch(ldSketch, jz$variantIds), + ldSketch = ldSketch, N = jz$nVec ) } @@ -1384,7 +1387,6 @@ setMethod("construct", "TwasJointPipeline", function(pipeline, records, ...) { ldSketch = ldSketch, cutoffs = args$cutoffs ) - ldMat <- .fmLdFromSketch(ldSketch, jz$variantIds) groups[[length(groups) + 1L]] <- new( "SumStatsJointGroup", conditions = tibble( @@ -1393,7 +1395,7 @@ setMethod("construct", "TwasJointPipeline", function(pipeline, records, ...) { trait = gi$traitCol[gIdx] ), Z = jz$Z, - R = ldMat, + ldSketch = ldSketch, N = jz$nVec ) } diff --git a/R/ld.R b/R/ld.R index 9eabb206..c2e54952 100644 --- a/R/ld.R +++ b/R/ld.R @@ -765,15 +765,17 @@ loadLdFromGenotype <- function( # ---------- LD sketch: per-variant LD matrix ---------- # Internal: build a sample-correlation LD matrix for a specified variant -# subset of an `ldSketch` `GenotypeHandle`. Shared by twasWeightsPipeline, -# fineMappingPipeline, causalInferencePipeline, and colocboostPipeline. The -# four sites differed only in their error message prefix and in whether -# variants absent from the panel raise an error or get silently dropped. +# subset of an `ldSketch` `GenotypeHandle`. Called directly by +# twasWeightsPipeline, fineMappingPipeline, jointEngine, +# causalInferencePipeline, colocboostPipeline and the coloc / summaryStatsQc +# paths, which differ only in their error message prefix and in whether +# variants absent from the panel raise an error or get silently dropped -- +# both of which are arguments here rather than a per-pipeline wrapper. # # Arguments: # ldSketch A GenotypeHandle. # variantIds Character vector of SNP IDs to extract. -# label Error-message prefix, e.g. ".twasLdFromSketch". +# label Error-message prefix, e.g. "twasWeightsPipeline". # onMissing "error" (default) -> any unmatched id stops the call; # "drop" -> unmatched ids are silently filtered. # When "drop" leaves no surviving variants the function @@ -890,10 +892,85 @@ loadLdFromGenotype <- function( .ldSketchHandle(x) } -# The panel's variant table, however the sketch is carried. +# ---------- LD sketch: the panel's own view of itself ---------- +# +# These read the panel through the SummarizedExperiment interface rather than +# reaching past it to the GenotypeHandle. The handle is the DelayedArray's +# seed -- the READ PATH -- and unwrapping it for metadata leaks that seed into +# callers that only ever wanted chrom/pos/alleles or a sample count, which the +# RSE already answers. Reaching for `.ldSketchHandle()` is now reserved for +# what genuinely lives on the handle: file format/path, chromosome shard +# routing, the `.afreq` sidecar, and building or pruning the seed itself. +# +# A bare GenotypeHandle -- what objects serialised before the slot became an +# RSE carry -- is projected exactly the way `.genotypeExperiment()` projects +# one, so both shapes answer identically. + +# The panel's variant ranges: seqnames = chromosome (chr-prefixed), start = +# position, mcols = SNP / A1 / A2. +# @noRd +.ldSketchRanges <- function(x) { + if (methods::is(x, "RangedSummarizedExperiment")) { + return(SummarizedExperiment::rowRanges(x)) + } + .genotypeSnpRanges(x, normalizeVariantId(as.character(getSnpInfo(x)$SNP))) +} + +# The panel's variant ids in the panel's OWN labelling. This is the mcols +# column, not `names()`: the names are normalized ids, while matching against +# summary statistics has to see the panel's raw labels. +# @noRd +.ldSketchVariantIds <- function(x) { + as.character(S4Vectors::mcols(.ldSketchRanges(x))$SNP) +} + +# Per-variant chromosome, canonical (no chr prefix). +# @noRd +.ldSketchChrom <- function(x) { + canonChrom(as.character(GenomicRanges::seqnames(.ldSketchRanges(x)))) +} + +# @noRd +.ldSketchSampleIds <- function(x) { + if (methods::is(x, "RangedSummarizedExperiment")) { + return(as.character(colnames(x))) + } + as.character(getSampleIds(x)) +} + +# @noRd +.ldSketchNSamples <- function(x) { + length(.ldSketchSampleIds(x)) +} + +# The BP / A1 / A2 vectors under the names the panel-comparison error messages +# report, so the message can name the column a user would recognise. +# @noRd +.ldRangesColumns <- function(gr) { + mc <- S4Vectors::mcols(gr) + list( + BP = as.character(GenomicRanges::start(gr)), + A1 = as.character(mc$A1), + A2 = as.character(mc$A2) + ) +} + +# Dosage for a variant subset, samples x variants. +# +# Read through the assay, so the DelayedArray seed -- and with it the file +# ordering, the sharded routing and the empty-request guard in +# `extractBlockGenotypes()` -- stays the single reader. `extract_array()` +# delegates to that same function with `meanImpute = FALSE`, so imputation is +# layered on here rather than pushed down. # @noRd -.ldSketchSnpInfo <- function(x) { - getSnpInfo(.ldSketchHandle(x)) +.ldSketchDosage <- function(x, snpIdx, meanImpute = TRUE) { + if (!methods::is(x, "RangedSummarizedExperiment")) { + return(.dosageMatrix(x, snpIdx, meanImpute = meanImpute)) + } + dosage <- t(as.matrix( + SummarizedExperiment::assay(x, "dosage")[snpIdx, , drop = FALSE] + )) + if (meanImpute) .qtlMeanImpute(dosage) else dosage } .ldFromSketchValidate <- function(ldSketch, label) { @@ -917,12 +994,11 @@ loadLdFromGenotype <- function( # Match requested ids to the panel; NULL if none match. Returns kept ids/order. .ldFromSketchMatch <- function(ldSketch, variantIds, label, onMissing) { - snpInfo <- .ldSketchSnpInfo(ldSketch) # Match by (chrom, pos, allele) tuple with an exact id-string fallback for # rsID panels; the caller's original ids and order are preserved. m <- matchVariants( variantIds, - as.character(snpInfo$SNP), + .ldSketchVariantIds(ldSketch), removeStrandAmbiguous = FALSE ) nMissing <- length(variantIds) - length(m$idxA) @@ -1060,7 +1136,7 @@ loadLdFromGenotype <- function( # would have to guess), so an imissCutoff always takes the dosage path. # @noRd .panelDropMask <- function( - handle, + ldSketch, matched, mafCutoff, macCutoff, @@ -1070,18 +1146,20 @@ loadLdFromGenotype <- function( # Keyed on the PANEL's own variant labels, not the caller's: .afreq # carries the .pvar ids, while `keptIds` is whatever id form the # caller asked in (normalized, most often). - panelIds <- as.character(getSnpInfo(handle)$SNP)[matched$idx] - maf <- .panelAfreqMaf(handle, panelIds) + panelIds <- .ldSketchVariantIds(ldSketch)[matched$idx] + # The handle stays for this one read: the sidecar is a file beside + # the genotypes, which only the seed knows how to find. + maf <- .panelAfreqMaf(.ldSketchHandle(ldSketch), panelIds) if (!is.null(maf)) { effMaf <- .panelEffectiveMaf( mafCutoff, macCutoff, - getNSamples(handle) + .ldSketchNSamples(ldSketch) ) return(is.na(maf) | maf < effMaf) } } - dosage <- .dosageMatrix(handle, matched$idx, meanImpute = FALSE) + dosage <- .ldSketchDosage(ldSketch, matched$idx, meanImpute = FALSE) stats <- .panelVariantStats(dosage) effMaf <- .panelEffectiveMaf(mafCutoff, macCutoff, nrow(dosage)) is.na(stats$maf) | stats$maf < effMaf | stats$missRate > imissCutoff @@ -1117,7 +1195,7 @@ loadLdFromGenotype <- function( return(variantIds) } drop <- .panelDropMask( - .ldSketchHandle(ldSketch), + ldSketch, matched, mafCutoff, macCutoff, @@ -1221,11 +1299,9 @@ loadLdFromGenotype <- function( ) abort(msg) } - qSnp <- .ldSketchSnpInfo(qtlLd) - gSnp <- .ldSketchSnpInfo(gwasLd) - if (nrow(qSnp) != nrow(gSnp)) { - nQ <- nrow(qSnp) - nG <- nrow(gSnp) + nQ <- length(.ldSketchRanges(qtlLd)) + nG <- length(.ldSketchRanges(gwasLd)) + if (nQ != nG) { msg <- glue( "{pipelineName}: ldSketch panels differ in size ({nQ} vs ", "{nG} variants){between}; the two ldSketch panels ", @@ -1237,17 +1313,19 @@ loadLdFromGenotype <- function( # Panels must agree on CHR/BP/A1/A2 columns and on the sample set. .ldSketchCheckContent <- function(qtlLd, gwasLd, pipelineName, between) { - qSnp <- .ldSketchSnpInfo(qtlLd) - gSnp <- .ldSketchSnpInfo(gwasLd) - if (!identical(canonChrom(qSnp$CHR), canonChrom(gSnp$CHR))) { + qGr <- .ldSketchRanges(qtlLd) + gGr <- .ldSketchRanges(gwasLd) + if (!identical(.ldSketchChrom(qtlLd), .ldSketchChrom(gwasLd))) { msg <- glue( "{pipelineName}: ldSketch panels differ in column CHR", "{between}; use the same ldSketch on both." ) abort(msg) } + qCol <- .ldRangesColumns(qGr) + gCol <- .ldRangesColumns(gGr) for (col in c("BP", "A1", "A2")) { - if (!identical(as.character(qSnp[[col]]), as.character(gSnp[[col]]))) { + if (!identical(qCol[[col]], gCol[[col]])) { msg <- glue( "{pipelineName}: ldSketch panels differ in column ", "{col}{between}; use the same ldSketch on both." @@ -1255,8 +1333,8 @@ loadLdFromGenotype <- function( abort(msg) } } - qIds <- getSampleIds(.ldSketchHandle(qtlLd)) - gIds <- getSampleIds(.ldSketchHandle(gwasLd)) + qIds <- .ldSketchSampleIds(qtlLd) + gIds <- .ldSketchSampleIds(gwasLd) if (!identical(qIds, gIds)) { msg <- glue( "{pipelineName}: ldSketch panels have different sample sets", @@ -2919,17 +2997,15 @@ loadLdBlock <- function(spec, g) { trimSamples, shrinkage ) { - handle <- .ldSketchHandle(panel) idx <- if (is.null(snpIdx)) { - seq_len(nrow(getSnpInfo(handle))) + seq_along(.ldSketchRanges(panel)) } else { snpIdx } - block <- extractBlockGenotypes(handle, idx) - if (is.null(block)) { + geno <- .ldSketchDosage(panel, idx) + if (is.null(geno)) { return(diag(length(idx))) } - geno <- t(SummarizedExperiment::assay(block, "dosage")) if (ncol(geno) < 2L) { return(diag(length(idx))) } @@ -2986,13 +3062,14 @@ loadLdBlock <- function(spec, g) { "'{fmt}'. Only GDS exposes an on-disk LD routine." )) } - handle <- .ldSketchHandle(X) idx <- if (is.null(snpIdx)) { - seq_len(nrow(getSnpInfo(handle))) + seq_along(.ldSketchRanges(X)) } else { snpIdx } - R <- .computeBlockLdGds(handle, idx) + # The handle stays here on purpose: the GDS on-disk LD routine reads the + # file directly, which is a seed-level operation with no assay equivalent. + R <- .computeBlockLdGds(.ldSketchHandle(X), idx) diag(R) <- 1.0 R[is.na(R) | is.nan(R)] <- 0 if (shrinkage > 0 && shrinkage <= 1) { diff --git a/R/manifestLoaders.R b/R/manifestLoaders.R index 85a69262..30719747 100644 --- a/R/manifestLoaders.R +++ b/R/manifestLoaders.R @@ -313,7 +313,7 @@ NULL if (length(getChromPaths(handle)) > 0L) { canonChrom(names(getChromPaths(handle))) } else { - unique(canonChrom(as.character(.ldSketchSnpInfo(ldSketch)$CHR))) + unique(.ldSketchChrom(ldSketch)) } } @@ -327,8 +327,8 @@ NULL return(invisible(NULL)) } reqChrom <- .ldContainmentCheckChroms(entryGr, ldSketch, label) - snpInfo <- .ldSketchSnpInfo(ldSketch) - if (is.null(snpInfo) || nrow(snpInfo) == 0L) { + sketchGr <- .ldSketchRanges(ldSketch) + if (length(sketchGr) == 0L) { msg <- glue( "{label}: LD sketch carries no variant metadata; ", "skipping the variant-overlap check." @@ -338,7 +338,7 @@ NULL } .ldContainmentCheckOverlap( entryGr, - snpInfo, + sketchGr, reqChrom, minLdOverlapWarn, label @@ -368,26 +368,21 @@ NULL # @noRd .ldContainmentCheckOverlap <- function( entryGr, - snpInfo, + sketchGr, reqChrom, minLdOverlapWarn, label ) { - keep <- is_in(canonChrom(as.character(snpInfo$CHR)), reqChrom) - snpInfo <- snpInfo[keep, , drop = FALSE] - entryVid <- formatVariantId( - chrom = canonChrom(as.character(GenomicRanges::seqnames(entryGr))), - pos = GenomicRanges::start(entryGr), - A2 = as.character(entryGr$A2), - A1 = as.character(entryGr$A1) + keep <- is_in( + canonChrom(as.character(GenomicRanges::seqnames(sketchGr))), + reqChrom ) - sketchVid <- formatVariantId( - chrom = canonChrom(as.character(snpInfo$CHR)), - pos = as.integer(snpInfo$BP), - A2 = as.character(snpInfo$A2), - A1 = as.character(snpInfo$A1) + nOverlap <- length( + matchVariants( + .grVariantIds(entryGr), + .grVariantIds(sketchGr[keep]) + )$idxA ) - nOverlap <- length(matchVariants(entryVid, sketchVid)$idxA) if (nOverlap == 0L) { msg <- glue( "{label}: none of the {length(entryGr)} summary-statistic ", diff --git a/R/overlapTopLoci.R b/R/overlapTopLoci.R index ad956495..2fbf8c27 100644 --- a/R/overlapTopLoci.R +++ b/R/overlapTopLoci.R @@ -162,7 +162,13 @@ setMethod( .overlapRelabelGwas <- function(gwasTl, vmap, coordCols) { g <- inner_join(gwasTl, vmap, by = c("variant_id" = "gwas_vid")) signedCols <- c("beta", "z", "conditional_effect") - g <- mutate(g, across(any_of(signedCols), .overlapApplySign, g$.sign)) + # Bound with partial() rather than across()'s deprecated `...`; the local + # matters because partial() resolves it lazily (see .negateWhere et al). + rowSign <- g$.sign + g <- mutate( + g, + across(any_of(signedCols), partial(.overlapApplySign, sign = rowSign)) + ) if (is_in("af", names(g))) { g$af <- if_else(g$.sign < 0 & !is.na(g$af), 1 - g$af, g$af) } diff --git a/R/qtlSumStats.R b/R/qtlSumStats.R index f24b2003..2dd4026e 100644 --- a/R/qtlSumStats.R +++ b/R/qtlSumStats.R @@ -538,8 +538,8 @@ setMethod("show", "QtlSumStats", function(object) { return(gr) } tp <- traitPos[i] - tssPos <- suppressWarnings(GenomicRanges::start(tp)) - tesPos <- suppressWarnings(GenomicRanges::end(tp)) + tssPos <- GenomicRanges::start(tp) + tesPos <- GenomicRanges::end(tp) if (length(tssPos) == 0L || is.na(tssPos)) { return(gr) } diff --git a/R/sumstatsQc.R b/R/sumstatsQc.R index c3d0c03f..56a6aca4 100644 --- a/R/sumstatsQc.R +++ b/R/sumstatsQc.R @@ -1727,18 +1727,37 @@ autoDecision <- function(df, highCorrCols) { } # Rows whose own allele flip is also present in the panel. +# +# "Flip" is the relation matchVariants() uses, not a literal A1/A2 reversal: +# a panel that records the second orientation on the opposite strand (A:G and +# C:T) carries exactly the ambiguity one that records it plainly (A:G and G:A) +# does, and a string key sees only the second. Since the mask exists to stop +# imputation putting back a variant matchVariants() dropped, the two have to +# mean the same thing -- so the panel is matched against ITSELF and a row is +# ambiguous precisely when a sumstats variant sitting on it would not survive. +# +# Self-matching also subsumes the degenerate A1 == A2 row (its own flip) that +# the string key had to special-case: such a row matches itself exactly and +# never as a swap, so no sign conflict arises. Unparseable ids fall back to +# exact string identity, where every row matches itself and nothing is +# ambiguous. +# +# The SECOND of two identical entries comes back TRUE, which the string key +# missed. That is not ambiguity but it is still unsafe: imputing it would put +# one id into the sumstats twice, and matchVariants answers a repeated id +# once, so the next LD lookup would abort on the copy it could not place. # @noRd .raissFlipPairMask <- function(refPanelIds) { n <- length(refPanelIds) - df <- parseVariantId(refPanelIds) - if (nrow(df) != n || anyNA(df$chrom) || anyNA(df$pos)) { - return(rep(FALSE, n)) + if (n == 0L) { + return(logical(0)) } - key <- str_c(df$chrom, df$pos, df$A2, df$A1, sep = ":") - flipped <- str_c(df$chrom, df$pos, df$A1, df$A2, sep = ":") - # A1 != A2 guard: were they equal a row would be its own flip and every - # variant would look ambiguous. - is_in(flipped, key) & df$A1 != df$A2 + m <- matchVariants( + refPanelIds, + refPanelIds, + removeStrandAmbiguous = FALSE + ) + !is_in(seq_len(n), m$idxA) } # Rows that match a known sumstats variant in the opposite orientation. @@ -2828,17 +2847,16 @@ krigingOutlierQc <- function( } # Build a refVariants data.frame (chrom, pos, A1, A2, variant_id) from the -# ldSketch GenotypeHandle's snpInfo so harmonizeAlleles can join by (chrom, -# pos). -.refVariantsFromSketch <- function(handle) { - si <- getSnpInfo(handle) - chr <- str_remove(as.character(si$CHR), regex("^chr", ignore_case = TRUE)) +# panel's own variant ranges so harmonizeAlleles can join by (chrom, pos). +.refVariantsFromSketch <- function(ldSketch) { + gr <- .ldSketchRanges(ldSketch) + mc <- S4Vectors::mcols(gr) data.frame( - chrom = chr, - pos = as.integer(si$BP), - A1 = as.character(si$A1), - A2 = as.character(si$A2), - variant_id = as.character(si$SNP), + chrom = .ldSketchChrom(ldSketch), + pos = as.integer(GenomicRanges::start(gr)), + A1 = as.character(mc$A1), + A2 = as.character(mc$A2), + variant_id = as.character(mc$SNP), stringsAsFactors = FALSE ) } @@ -3141,7 +3159,7 @@ krigingOutlierQc <- function( removeStrandAmbiguous = TRUE, removeDups = TRUE ) { - refVariants <- .refVariantsFromSketch(.ldSketchHandle(ldSketch)) + refVariants <- .refVariantsFromSketch(ldSketch) flipCandidates <- c("Z", "BETA") colToFlip <- intersect(flipCandidates, colnames(df)) if (length(colToFlip) == 0L) { @@ -3529,7 +3547,7 @@ krigingOutlierQc <- function( qc <- ldMismatchQc( zScore = df$Z, R = R, - nSample = getNSamples(.ldSketchHandle(ldSketch)), + nSample = .ldSketchNSamples(ldSketch), method = method ) # slalom / dentist can leave NA in the outlier column when their @@ -3817,12 +3835,11 @@ krigingOutlierQc <- function( # --- .runEntrySummaryStatsQc: RAISS imputation step helpers ---------------- # Panel/dosage window indices for the entry, scoped per chromosome. -.qcRaissWindowIdx <- function(df, sketchSnpInfo, flank) { +.qcRaissWindowIdx <- function(df, ldSketch, flank) { + # canonChrom on BOTH sides: the panel answers in its own seqnames + # convention, so the two are only comparable once canonicalized. bounds <- tibble( - chrom = str_remove( - as.character(df$chrom), - regex("^chr", ignore_case = TRUE) - ), + chrom = canonChrom(as.character(df$chrom)), pos = as.integer(df$pos) ) |> group_by(.data$chrom) |> @@ -3834,13 +3851,11 @@ krigingOutlierQc <- function( # inner_join keeps only sketch SNPs whose chromosome appears in df (the # old is_in(skChrom, names(loByChr)) guard); filter keeps those inside the # [lo, hi] window. idx carries the original sketch row positions. + gr <- .ldSketchRanges(ldSketch) tibble( - chrom = str_remove( - as.character(sketchSnpInfo$CHR), - regex("^chr", ignore_case = TRUE) - ), - bp = as.integer(sketchSnpInfo$BP), - idx = seq_along(sketchSnpInfo$CHR) + chrom = .ldSketchChrom(ldSketch), + bp = as.integer(GenomicRanges::start(gr)), + idx = seq_along(gr) ) |> inner_join(bounds, by = "chrom") |> filter(.data$bp >= .data$lo & .data$bp <= .data$hi) |> @@ -3879,23 +3894,17 @@ krigingOutlierQc <- function( arrange(knownZ, .data$pos) } -.qcRaissBuildInputs <- function(df, ldSketch, windowIdx, sketchSnpInfo, opts) { - refPanel <- .refVariantsFromSketch( - .ldSketchHandle(ldSketch) - )[windowIdx, , drop = FALSE] +.qcRaissBuildInputs <- function(df, ldSketch, windowIdx, opts) { + refPanel <- .refVariantsFromSketch(ldSketch)[windowIdx, , drop = FALSE] refPanel$variant_id <- normalizeVariantId(refPanel$variant_id) refPanel <- arrange(refPanel, .data$pos) knownZ <- .qcRaissKnownZ(df) # meanImpute = FALSE so per-variant missingness is still visible; the # surviving columns are mean-imputed below, which is what meanImpute = # TRUE did. - dosage <- .dosageMatrix( - .ldSketchHandle(ldSketch), - windowIdx, - meanImpute = FALSE - ) + dosage <- .ldSketchDosage(ldSketch, windowIdx, meanImpute = FALSE) colnames(dosage) <- normalizeVariantId( - as.character(sketchSnpInfo$SNP[windowIdx]) + .ldSketchVariantIds(ldSketch)[windowIdx] ) dosage <- dosage[, refPanel$variant_id, drop = FALSE] keep <- .qcRaissTargetMask(refPanel, knownZ, dosage, opts) @@ -4054,8 +4063,7 @@ krigingOutlierQc <- function( } else { as.integer(opts$imputeOpts$flank) } - sketchSnpInfo <- .ldSketchSnpInfo(ldSketch) - windowIdx <- .qcRaissWindowIdx(df, sketchSnpInfo, flank) + windowIdx <- .qcRaissWindowIdx(df, ldSketch, flank) if (length(windowIdx) == 0L) { .qcEmit( lbl, @@ -4071,7 +4079,7 @@ krigingOutlierQc <- function( imputeAfter = nrow(df) )) } - inp <- .qcRaissBuildInputs(df, ldSketch, windowIdx, sketchSnpInfo, opts) + inp <- .qcRaissBuildInputs(df, ldSketch, windowIdx, opts) if (inp$nDroppedTargets > 0L) { .qcEmit( lbl, @@ -4557,14 +4565,14 @@ krigingOutlierQc <- function( bounds <- tibble(chrom = chrom, pos = pos) |> group_by(chrom) |> summarise(lo = min(pos), hi = max(pos), .groups = "drop") - si <- .ldSketchSnpInfo(ldSketch) + gr <- .ldSketchRanges(ldSketch) # left_join keeps every sketch SNP; one on a chromosome absent from the # entries gets lo/hi = NA -> inWindow FALSE (the old is_in guard). The # explicit is.na() guards force FALSE (never NA) for an absent chrom or an # NA bp, replacing the base keep[is.na(keep)] <- FALSE. keep <- tibble( - chrom = canonChrom(as.character(si$CHR)), - bp = as.integer(si$BP) + chrom = .ldSketchChrom(ldSketch), + bp = as.integer(GenomicRanges::start(gr)) ) |> left_join(bounds, by = "chrom") |> mutate( @@ -4599,9 +4607,8 @@ krigingOutlierQc <- function( if (length(ids) == 0L) { return(ldSketch) } - si <- .ldSketchSnpInfo(ldSketch) keep <- is_in( - normalizeVariantId(as.character(si$SNP)), + normalizeVariantId(.ldSketchVariantIds(ldSketch)), normalizeVariantId(unique(ids)) ) .ldSketchSubset(ldSketch, keep) @@ -4751,8 +4758,7 @@ krigingOutlierQc <- function( if (is.null(ldSketch) || is.null(cutoffs)) { return(ldSketch) } - handle <- .ldSketchHandle(ldSketch) - ids <- as.character(getSnpInfo(handle)$SNP) + ids <- .ldSketchVariantIds(ldSketch) if (length(ids) == 0L) { return(ldSketch) } @@ -4760,12 +4766,12 @@ krigingOutlierQc <- function( if (all(keep)) { return(ldSketch) } - # Prune the underlying GenotypeHandle, not just an RSE row view. Every - # reader reaches the panel through `.ldSketchHandle()` -- the dosage - # DelayedArray's seed -- and `sketch[keep, ]` leaves that seed carrying - # the whole panel (the trap `.emptySketch()` documents). Rebuild the - # sketch from the pruned handle so the drop reaches all of them. - pruned <- .subsetGenotypeHandle(handle, keep) + # Prune the underlying GenotypeHandle, not just an RSE row view: the + # handle is the dosage DelayedArray's seed, and `sketch[keep, ]` leaves + # that seed carrying the whole panel (the trap `.emptySketch()` + # documents), so every read would still see the dropped variants. This is + # a seed-level edit, which is why it reaches for the handle. + pruned <- .subsetGenotypeHandle(.ldSketchHandle(ldSketch), keep) if (methods::is(ldSketch, "GenotypeHandle")) { return(pruned) } diff --git a/R/twasWeightsPipeline.R b/R/twasWeightsPipeline.R index 1991d65a..3a185e53 100644 --- a/R/twasWeightsPipeline.R +++ b/R/twasWeightsPipeline.R @@ -743,14 +743,6 @@ combineTwasWeights <- function(..., ldSketch = NULL) { } } -# Extract a correlation matrix from a GenotypeHandle (LD sketch) for the -# variant subset given by `variantIds`. Thin wrapper over the shared -# `.ldFromSketch` helper. -# @noRd -.twasLdFromSketch <- function(ldSketch, variantIds) { - .ldFromSketch(ldSketch, variantIds, label = ".twasLdFromSketch") -} - # Optional resume-cache lookup for twasWeightsPipeline. Returns the # matching TwasWeightsRow from `twasWeights` for the tuple (study, # context, trait, method), or NULL when there is no hit. Returns NULL @@ -1862,7 +1854,11 @@ setMethod( list( variantIds = variantIds, stat = stat, - ldMat = .twasLdFromSketch(ldSketch, variantIds) + ldMat = .ldFromSketch( + ldSketch, + variantIds, + label = "twasWeightsPipeline" + ) ) } @@ -1968,7 +1964,11 @@ setMethod( ldSketch = p$ldSketch, cutoffs = .panelCutoffs(p) ) - ldMat <- .twasLdFromSketch(p$ldSketch, mvStat$variantIds) + ldMat <- .ldFromSketch( + p$ldSketch, + mvStat$variantIds, + label = "twasWeightsPipeline" + ) list_flatten(map( p$multivariateTokens, .twasQssMultivariateFitOne, diff --git a/R/variantId.R b/R/variantId.R index 2177dfcf..b343a4ed 100644 --- a/R/variantId.R +++ b/R/variantId.R @@ -412,6 +412,7 @@ variantIdToDf <- function(variantId) { #' per-variant merge/flip/strand diagnostics. #' @importFrom dplyr mutate inner_join filter pull select everything row_number #' @importFrom dplyr if_else any_of all_of rename across +#' @importFrom purrr partial #' @importFrom vctrs vec_duplicate_detect #' @importFrom tidyr separate #' @keywords internal @@ -439,21 +440,20 @@ harmonizeAlleles <- function( coerced <- .harmonizeCoerceInputs(targetData, refVariants) targetData <- coerced$targetData refVariants <- coerced$refVariants - matchResult <- .harmonizeJoin(targetData, refVariants) + # The index rides along on the join copy only: the restore path below + # re-reads `targetData`, and a sentinel column there would bind_rows its + # way into the returned frame. + matchResult <- .harmonizeJoin(.haIndexTargets(targetData), refVariants) if (nrow(matchResult) == 0) { return(.harmonizeEmptyResult(matchResult)) } - matchResult <- .harmonizeFlags(matchResult) - matchResult <- .harmonizeResolveAmbiguity( - matchResult, - removeStrandAmbiguous - ) - matchResult <- .harmonizeKeepRule(matchResult, removeIndels) - matchResult <- .harmonizeApplyFlips( + matchResult <- .harmonizeDecideMatches( matchResult, colToFlip, colToComplement, - flipStrand + flipStrand, + removeIndels, + removeStrandAmbiguous ) qcCounts <- .harmonizeQcCounts(matchResult) qcSummary <- matchResult @@ -472,6 +472,48 @@ harmonizeAlleles <- function( out } +# Turn the raw (chrom, pos) join into a decided match table: classify each +# pair's allele relationship, drop the strand-ambiguous and indel rows the +# caller refused, settle each target on a single reference row, then apply the +# sign / frequency flips the chosen orientation implies. Everything before this +# builds the join; everything after it shapes the output. +# @noRd +.harmonizeDecideMatches <- function( + matchResult, + colToFlip, + colToComplement, + flipStrand, + removeIndels, + removeStrandAmbiguous +) { + matchResult <- .harmonizeFlags(matchResult) + matchResult <- .harmonizeResolveAmbiguity( + matchResult, + removeStrandAmbiguous + ) + matchResult <- .harmonizeKeepRule(matchResult, removeIndels) + matchResult <- .harmonizeResolveTargets(matchResult) + .harmonizeApplyFlips( + matchResult, + colToFlip, + colToComplement, + flipStrand + ) +} + +# Target-row index carried through the (chrom, pos) join so the resolution +# below can group the reference rows one target variant matched. +.haTargetIdx <- ".haTidx" + +# Assigned directly rather than through mutate(): the target frame is not +# name-repaired until the join (.sanitizeNames), and dplyr's data mask refuses +# a frame carrying an NA or empty column name. +# @noRd +.haIndexTargets <- function(targetData) { + targetData[[.haTargetIdx]] <- seq_len(nrow(targetData)) + targetData +} + # QC / flag columns stripped from the harmonized result before it is returned. .harmonizeQcCols <- c( "flip1.ref", @@ -482,7 +524,8 @@ harmonizeAlleles <- function( "strand_flip", "INDEL", "ID_match", - "keep" + "keep", + .haTargetIdx ) # Coerce both sides to canonical variant data.frames and strip merge-conflicting @@ -667,9 +710,66 @@ harmonizeAlleles <- function( matchResult } -# Named per-row conditional column transforms for the harmonize across() calls -# (the `flip` condition vector is passed through across's `...`, so no anonymous -# functions are needed). +# Reduce the join to at most one reference row per TARGET variant. +# +# The join is on (chrom, pos), so a target matches every reference row at its +# position whose alleles reconcile. Usually that is one row. It is more when +# the reference carries the variant AND its own flip as separate entries -- +# two distinct indels at one position that happen to be each other's flip is +# rare but biologically real -- or the same variant once per strand. +# +# Matches that DISAGREE about the allele-swap sign are undecidable: the target +# is either the exact entry or the swapped one and nothing says which. Keeping +# the first would make the effect direction depend on the reference's row +# order; keeping both would return the target twice with opposite signs, which +# is not harmonization but duplication. Those targets are dropped. Matches that +# agree (the same variant listed twice, or once per strand) resolve to the +# first row, so this costs nothing where there is no ambiguity. +# +# Deciding it here, before the keep flags are counted, is what lets every +# caller -- matchVariants, the summaryStatsQc panel harmonization, the cTWAS +# weight harmonization -- agree on the surviving variant set without repeating +# the rule. +# @noRd +.harmonizeResolveTargets <- function(matchResult) { + tidx <- matchResult[[.haTargetIdx]] + hit <- which(matchResult$keep) + # Fast path: one reference row per target is the overwhelmingly common + # case, and it needs no grouping at all. + if (anyDuplicated(tidx[hit]) == 0L) { + return(matchResult) + } + conflicted <- .harmonizeConflictedTargets( + tidx[hit], + matchResult$sign_flip[hit] + ) + resolved <- hit[ + !is_in(tidx[hit], conflicted) & !duplicated(tidx[hit]) + ] + # Named apart from the column it replaces: `keep = keep` inside mutate() + # reads the column, not this vector, and silently changes nothing. + resolvedKeep <- rep(FALSE, nrow(matchResult)) + resolvedKeep[resolved] <- TRUE + mutate(matchResult, keep = resolvedKeep) +} + +# Target indices whose reference matches disagree about the allele-swap sign. +# @noRd +.harmonizeConflictedTargets <- function(tidx, signFlip) { + tibble(tidx = tidx, signFlip = signFlip) |> + group_by(.data$tidx) |> + summarise(nSigns = n_distinct(.data$signFlip), .groups = "drop") |> + filter(.data$nSigns > 1L) |> + pull("tidx") +} + +# Named per-row conditional column transforms for the harmonize across() calls. +# The `flip` condition vector is bound with purrr::partial() rather than passed +# through across()'s `...`, which dplyr deprecated in 1.1.0 -- partial() fixes +# the argument without an anonymous function, so the no-lambda rule the `...` +# form was chosen for still holds. partial() is LAZY, so callers bind the +# condition to a local first: what it captures must not be reassigned before +# across() calls it. # @noRd .negateWhere <- function(x, flip) if_else(flip, -x, x) # @noRd @@ -686,13 +786,13 @@ harmonizeAlleles <- function( colToComplement, flipStrand ) { + signFlip <- matchResult$sign_flip if (!is.null(colToFlip)) { .harmonizeCheckCols(colToFlip, matchResult) matchResult <- matchResult |> mutate(across( all_of(colToFlip), - .negateWhere, - matchResult$sign_flip + partial(.negateWhere, flip = signFlip) )) } if (length(colToComplement) > 0L) { @@ -700,8 +800,7 @@ harmonizeAlleles <- function( matchResult <- matchResult |> mutate(across( all_of(colToComplement), - .complementWhere, - matchResult$sign_flip + partial(.complementWhere, flip = signFlip) )) } if (flipStrand) { @@ -725,11 +824,11 @@ harmonizeAlleles <- function( # Strand-flip the target alleles of the strand-flipped rows. # @noRd .harmonizeFlipStrandCols <- function(matchResult) { + strandFlip <- matchResult$strand_flip matchResult |> mutate(across( c("A1.target", "A2.target"), - .strandFlipWhere, - matchResult$strand_flip + partial(.strandFlipWhere, flip = strandFlip) )) } @@ -993,44 +1092,16 @@ matchVariants <- function( if (is.null(h) || nrow(h) == 0L) { return(.matchVariantsEmpty()) } - h <- .matchDropSignConflicts(h) - if (nrow(h) == 0L) { - return(.matchVariantsEmpty()) - } - keep <- !duplicated(h$.mvTidx) # at most one ref per A id + # One row per A id, sign-conflicting ids already dropped: harmonizeAlleles + # resolves that (see .harmonizeResolveTargets), so there is nothing to + # decide a second time here. list( - idxA = as.integer(h$.mvTidx[keep]), - idxB = as.integer(h$.mvRidx[keep]), - sign = as.numeric(h$.mvSign[keep]) + idxA = as.integer(h$.mvTidx), + idxB = as.integer(h$.mvRidx), + sign = as.numeric(h$.mvSign) ) } -# Drop target variants whose reference matches disagree about the sign. -# -# A reference panel can legitimately carry a variant and its own allele flip as -# two separate entries -- two distinct indels at one position that happen to be -# each other's flip is rare but biologically real. A target variant then -# matches BOTH: once exactly (sign +1) and once as an allele swap (sign -1), -# and there is no way to tell which entry it is. The caller below keeps the -# first match, so without this the answer would depend on the panel's row -# order, silently flipping the effect direction for that variant. -# -# Only a sign DISAGREEMENT is ambiguous. Duplicate reference entries that agree -# (the same variant listed twice) still resolve to the first match, so this -# drops nothing it does not have to. -# @noRd -.matchDropSignConflicts <- function(h) { - ambiguous <- h |> - group_by(.data$.mvTidx) |> - summarise(nSigns = n_distinct(.data$.mvSign), .groups = "drop") |> - filter(.data$nSigns > 1L) |> - pull(".mvTidx") - if (length(ambiguous) == 0L) { - return(h) - } - filter(h, !is_in(.data$.mvTidx, ambiguous)) -} - # Backwards-compat alias for external callers #' Parse a region string into its components diff --git a/man/LdData-class.Rd b/man/LdData-class.Rd index adcb75ad..28fd9614 100644 --- a/man/LdData-class.Rd +++ b/man/LdData-class.Rd @@ -3,12 +3,43 @@ \docType{class} \name{LdData-class} \alias{LdData-class} +\alias{[,LdData,ANY,ANY,ANY-method} \title{LD Data Container} +\usage{ +\S4method{[}{LdData,ANY,ANY,ANY}(x, i, j, ..., drop = TRUE) +} +\arguments{ +\item{x}{An \code{LdData}.} + +\item{i, j, ...}{Subscripts; any use is an error.} + +\item{drop}{Ignored.} +} +\value{ +Nothing: this method always signals an error. +} \description{ S4 container for LD information. Stores either a pre-computed correlation matrix or a \code{GenotypeHandle} (or list of handles for mixture panels) for lazy genotype/correlation access. + + An \code{LdData} \strong{is} a \code{GRanges} over the variants it + covers -- carrying A1, A2, variant_id and optionally allele_freq, + variance and n_nomiss as metadata columns -- so \code{length()}, + \code{seqnames()} and \code{start()} answer directly, as they do for + \code{\link{LdScore}} and \code{\link{LdEigen}}. The correlation may + be NULL, in which case those ranges are the object's only record of which + variants it describes. } +\section{Functions}{ +\itemize{ +\item \code{x = LdData[i = ANY, j = ANY, drop = ANY]}: Refused. Subsetting would narrow the variants +while \code{correlation} -- variant-by-variant, and \code{snpIdx}, which +indexes the reference panel -- stayed as they were, leaving an LD matrix +describing variants the object no longer has. Build the LdData over the +variant set you want instead. + +}} \section{Slots}{ \describe{ @@ -26,9 +57,6 @@ it is unwrapped to its handle.} \code{snpInfo}. NULL when correlation is pre-computed, or when the source is a matrix (which is already the subset).} -\item{\code{variants}}{A \code{GRanges} object with variant metadata (A1, A2, -variant_id, and optionally allele_freq, variance, n_nomiss).} - \item{\code{blockMetadata}}{A \code{GRanges} of blocks or a \code{data.frame} with block boundary information.} diff --git a/man/getH2Se.Rd b/man/getH2Se.Rd new file mode 100644 index 00000000..9cd9d1c4 --- /dev/null +++ b/man/getH2Se.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/AllGenerics.R, R/H2Estimate.R +\name{getH2Se} +\alias{getH2Se} +\alias{getH2Se,H2Estimate-method} +\title{Get Heritability Standard Error} +\usage{ +getH2Se(x) + +\S4method{getH2Se}{H2Estimate}(x) +} +\arguments{ +\item{x}{An \code{H2Estimate}.} +} +\value{ +Numeric (length 1). +} +\description{ +Return the standard error of the global SNP heritability + estimate carried by an \code{H2Estimate}. +} +\examples{ +data(h2EstimateExample) +getH2Se(h2EstimateExample) +} diff --git a/man/getIntercept.Rd b/man/getIntercept.Rd new file mode 100644 index 00000000..e037f2ff --- /dev/null +++ b/man/getIntercept.Rd @@ -0,0 +1,26 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/AllGenerics.R, R/H2Estimate.R +\name{getIntercept} +\alias{getIntercept} +\alias{getIntercept,H2Estimate-method} +\title{Get LD-Score Regression Intercept} +\usage{ +getIntercept(x) + +\S4method{getIntercept}{H2Estimate}(x) +} +\arguments{ +\item{x}{An \code{H2Estimate}.} +} +\value{ +Numeric (length 1). +} +\description{ +Return the regression intercept carried by an + \code{H2Estimate}. Values above 1 indicate confounding or sample overlap + rather than polygenic signal. +} +\examples{ +data(h2EstimateExample) +getIntercept(h2EstimateExample) +} diff --git a/man/getInterceptSe.Rd b/man/getInterceptSe.Rd new file mode 100644 index 00000000..0a84b72c --- /dev/null +++ b/man/getInterceptSe.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/AllGenerics.R, R/H2Estimate.R +\name{getInterceptSe} +\alias{getInterceptSe} +\alias{getInterceptSe,H2Estimate-method} +\title{Get Intercept Standard Error} +\usage{ +getInterceptSe(x) + +\S4method{getInterceptSe}{H2Estimate}(x) +} +\arguments{ +\item{x}{An \code{H2Estimate}.} +} +\value{ +Numeric (length 1). +} +\description{ +Return the standard error of the regression intercept carried + by an \code{H2Estimate}. +} +\examples{ +data(h2EstimateExample) +getInterceptSe(h2EstimateExample) +} diff --git a/man/getMethodNames.Rd b/man/getMethodNames.Rd index 2268cfcb..3cf5f3a5 100644 --- a/man/getMethodNames.Rd +++ b/man/getMethodNames.Rd @@ -1,10 +1,11 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/AllGenerics.R, R/AllClasses.R, -% R/CtwasResult.R, R/twasWeights.R +% R/CtwasResult.R, R/H2Estimate.R, R/twasWeights.R \name{getMethodNames} \alias{getMethodNames} \alias{getMethodNames,FineMappingResultBase-method} \alias{getMethodNames,CtwasResult-method} +\alias{getMethodNames,H2Estimate-method} \alias{getMethodNames,TwasWeights-method} \title{Get Method Names} \usage{ @@ -14,18 +15,24 @@ getMethodNames(x) \S4method{getMethodNames}{CtwasResult}(x) +\S4method{getMethodNames}{H2Estimate}(x) + \S4method{getMethodNames}{TwasWeights}(x) } \arguments{ -\item{x}{A \code{FineMappingResult} or \code{TwasWeights} object.} +\item{x}{A \code{FineMappingResult}, \code{TwasWeights}, +\code{CtwasResult} or \code{H2Estimate} object.} } \value{ Character vector. } \description{ -Extract method names from a collection class. +Extract the method name(s) an object was produced with: one + per entry for a collection, a single name for one estimate. } \examples{ +data(h2EstimateExample) +getMethodNames(h2EstimateExample) data(qtlFineMappingExample) getMethodNames(qtlFineMappingExample) } diff --git a/man/getNSnps.Rd b/man/getNSnps.Rd new file mode 100644 index 00000000..9f8d948f --- /dev/null +++ b/man/getNSnps.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/AllGenerics.R, R/H2Estimate.R +\name{getNSnps} +\alias{getNSnps} +\alias{getNSnps,H2Estimate-method} +\title{Get Variant Count} +\usage{ +getNSnps(x) + +\S4method{getNSnps}{H2Estimate}(x) +} +\arguments{ +\item{x}{An \code{H2Estimate}.} +} +\value{ +Integer (length 1). +} +\description{ +Return the number of variants the heritability estimate was + computed over. +} +\examples{ +data(h2EstimateExample) +getNSnps(h2EstimateExample) +} diff --git a/man/getSnpInfo.Rd b/man/getSnpInfo.Rd index 5807d270..b2837dc9 100644 --- a/man/getSnpInfo.Rd +++ b/man/getSnpInfo.Rd @@ -17,6 +17,15 @@ A data.frame. } \description{ Return the cached SNP metadata data.frame (columns: SNP, CHR, - BP, A1, A2, optionally MAF). + BP, A1, A2, fileIdx, optionally MAF). +} +\details{ +The genotype handle is the seed layer behind a panel and + is not part of the public interface; obtain a panel from + \code{readGenotypes()} and ask it directly. This table is the READ-PATH + view: unlike the panel's \code{rowRanges()} it carries \code{fileIdx}, + the on-disk variant position that keeps reads correct after a handle has + been row-subset. Callers that only need chrom/pos/alleles should use + \code{rowRanges()} on the panel instead. } \keyword{internal} diff --git a/man/getTraitName.Rd b/man/getTraitName.Rd new file mode 100644 index 00000000..f6607d2c --- /dev/null +++ b/man/getTraitName.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/AllGenerics.R, R/H2Estimate.R +\name{getTraitName} +\alias{getTraitName} +\alias{getTraitName,H2Estimate-method} +\title{Get Trait Name} +\usage{ +getTraitName(x) + +\S4method{getTraitName}{H2Estimate}(x) +} +\arguments{ +\item{x}{An \code{H2Estimate}.} +} +\value{ +Character (length 1). +} +\description{ +Return the trait label carried by an \code{H2Estimate}. +} +\examples{ +data(h2EstimateExample) +getTraitName(h2EstimateExample) +} diff --git a/tests/testthat/helper-h2Classes.R b/tests/testthat/helper-h2Classes.R index e261e5ac..9d3418ac 100644 --- a/tests/testthat/helper-h2Classes.R +++ b/tests/testthat/helper-h2Classes.R @@ -143,8 +143,13 @@ makeGwasSumStatsFromDf <- function( # One LD block's worth of GWAS summary statistics: `n` variants starting at # `start`, carrying a QC record and an LD panel narrowed to just those # variants (what a per-block pipeline step writes out). -makeGwasBlock <- function(blockId, start, n = 5L, genome = "hg19", - qcOptions = list(mafCutoff = 0.01)) { +makeGwasBlock <- function( + blockId, + start, + n = 5L, + genome = "hg19", + qcOptions = list(mafCutoff = 0.01) +) { bp <- seq(start, by = 100L, length.out = n) df <- data.frame( SNP = paste0("chr1:", bp, ":A:G"), diff --git a/tests/testthat/test_ColocBoostResult.R b/tests/testthat/test_ColocBoostResult.R index e6de3a0d..8aa18856 100644 --- a/tests/testthat/test_ColocBoostResult.R +++ b/tests/testthat/test_ColocBoostResult.R @@ -304,7 +304,9 @@ test_that("show on an empty ColocBoostResult stops after the header", { x <- ColocBoostResult(list(), character(0), outcomeInfo = .cbr_info()) expect_output(show(x), "with 0 confidence set\\(s\\)") expect_false(any(grepl( - "max cos_npc", capture.output(show(x)), fixed = TRUE + "max cos_npc", + capture.output(show(x)), + fixed = TRUE ))) }) diff --git a/tests/testthat/test_ColocResult.R b/tests/testthat/test_ColocResult.R index d33d1f37..215021fd 100644 --- a/tests/testthat/test_ColocResult.R +++ b/tests/testthat/test_ColocResult.R @@ -254,7 +254,9 @@ test_that("show on an empty ColocResult stops after the header", { x <- ColocResult(.cr_pairs()[0, , drop = FALSE], list()) expect_output(show(x), "with 0 colocalized pair\\(s\\)") expect_false(any(grepl( - "max PP.H4", capture.output(show(x)), fixed = TRUE + "max PP.H4", + capture.output(show(x)), + fixed = TRUE ))) }) diff --git a/tests/testthat/test_GwasSumStats.R b/tests/testthat/test_GwasSumStats.R index 9686179d..7634bd6e 100644 --- a/tests/testthat/test_GwasSumStats.R +++ b/tests/testthat/test_GwasSumStats.R @@ -450,12 +450,17 @@ test_that("GwasSumStats blockId survives subsetting", { .gss_entry <- function(n = 3L) { gr <- GenomicRanges::GRanges( seqnames = rep("chr1", n), - ranges = IRanges::IRanges(seq(100L, by = 100L, length.out = n), - width = 1L) + ranges = IRanges::IRanges( + seq(100L, by = 100L, length.out = n), + width = 1L + ) ) S4Vectors::mcols(gr) <- S4Vectors::DataFrame( - variant_id = str_c("chr1:", seq(100L, by = 100L, length.out = n), - ":A:G"), + variant_id = str_c( + "chr1:", + seq(100L, by = 100L, length.out = n), + ":A:G" + ), SNP = str_c("rs", seq_len(n)), A1 = rep("A", n), A2 = rep("G", n), @@ -472,7 +477,8 @@ test_that("supplying both ldBlocks and blockId is refused", { entry = list(.gss_entry()), genome = "hg19", ldBlocks = GenomicRanges::GRanges( - "chr1", IRanges::IRanges(1L, 10000L) + "chr1", + IRanges::IRanges(1L, 10000L) ), blockId = "b1" ), @@ -570,13 +576,21 @@ test_that("combineGwasSumStats() rejects collection-level disagreement", { a <- makeGwasBlock("chr1_1_1000", 100L) expect_error( - combineGwasSumStats(a, makeGwasBlock("chr1_1001_2000", 1100L, - genome = "hg38")), + combineGwasSumStats( + a, + makeGwasBlock("chr1_1001_2000", 1100L, genome = "hg38") + ), "share one genome build" ) expect_error( - combineGwasSumStats(a, makeGwasBlock("chr1_1001_2000", 1100L, - qcOptions = list(mafCutoff = 0.05))), + combineGwasSumStats( + a, + makeGwasBlock( + "chr1_1001_2000", + 1100L, + qcOptions = list(mafCutoff = 0.05) + ) + ), "different summaryStatsQc\\(\\) options" ) noQc <- makeGwasBlock("chr1_1001_2000", 1100L) @@ -589,8 +603,10 @@ test_that("combineGwasSumStats() rejects collection-level disagreement", { other <- makeGwasBlock("chr1_1001_2000", 1100L) other@ldSketch <- pecotmr:::.asLdSketch( - .blockGenotypeHandle(seq(1100L, by = 100L, length.out = 5L), - path = "/tmp/other.pgen") + .blockGenotypeHandle( + seq(1100L, by = 100L, length.out = 5L), + path = "/tmp/other.pgen" + ) ) expect_error( combineGwasSumStats(a, other), diff --git a/tests/testthat/test_H2Estimate.R b/tests/testthat/test_H2Estimate.R index 1e139652..be452809 100644 --- a/tests/testthat/test_H2Estimate.R +++ b/tests/testthat/test_H2Estimate.R @@ -18,8 +18,8 @@ test_that("H2Estimate constructs with all slots", { traitName = "height" ) expect_s4_class(obj, "H2Estimate") - expect_equal(obj@h2, 0.3) - expect_equal(obj@method, "lder") + expect_equal(getH2(obj), 0.3) + expect_equal(getMethodNames(obj), "lder") }) # show() smoke test, moved here from test_showMethods.R so the test diff --git a/tests/testthat/test_JointGroup.R b/tests/testthat/test_JointGroup.R index 03da66d2..4531dc38 100644 --- a/tests/testthat/test_JointGroup.R +++ b/tests/testthat/test_JointGroup.R @@ -23,15 +23,14 @@ test_that("JointGroup subclasses construct from a conditions table", { g <- new("IndividualJointGroup", conditions = .jg_cond(), X = X, Y = Y) expect_s4_class(g, "JointGroup") expect_s4_class(g, "IndividualJointGroup") - expect_equal(nrow(g@conditions), 2L) + expect_equal(nrow(.jgConditions(g)), 2L) - Z <- matrix(0, 3, 2) - R <- diag(3) + Z <- matrix(0, 3, 2, dimnames = list(paste0("v", 1:3), NULL)) sg <- new( "SumStatsJointGroup", conditions = .jg_cond(), Z = Z, - R = R, + ldSketch = NULL, N = c(100, 120) ) expect_s4_class(sg, "JointGroup") @@ -99,16 +98,17 @@ test_that("JointGroup validity rejects malformed groups", { ), "ncol\\(Y\\)" ) - # Non-square LD. + # Z must name its variants: they are what the LD matrix is derived over, + # so an unnamed Z would leave the group unable to say what it covers. expect_error( new( "SumStatsJointGroup", conditions = .jg_cond(), Z = matrix(0, 3, 2), - R = matrix(0, 3, 2), + ldSketch = NULL, N = 1 ), - "square" + "variant ids as rownames" ) }) diff --git a/tests/testthat/test_QtlDataset.R b/tests/testthat/test_QtlDataset.R index e1017abc..8568e4f8 100644 --- a/tests/testthat/test_QtlDataset.R +++ b/tests/testthat/test_QtlDataset.R @@ -5,7 +5,14 @@ context("QtlDataset internal helpers") # functions are stubbed in individual tests via local_mocked_bindings. # =========================================================================== -.qh_makeHandle <- function(snp_n = 6L, n_samples = 12L) { +# `a1` / `a2` are parameters so a test needing indels asks for them, rather +# than reaching into the built handle's snpInfo to rewrite the alleles. +.qh_makeHandle <- function( + snp_n = 6L, + n_samples = 12L, + a1 = rep("A", snp_n), + a2 = rep("G", snp_n) +) { new( "GenotypeHandle", path = "/tmp/test.gds", @@ -14,8 +21,8 @@ context("QtlDataset internal helpers") SNP = paste0("rs", seq_len(snp_n)), CHR = rep("1", snp_n), BP = seq(100L, by = 100L, length.out = snp_n), - A1 = rep("A", snp_n), - A2 = rep("G", snp_n), + A1 = a1, + A2 = a2, stringsAsFactors = FALSE ), nSamples = n_samples, @@ -63,12 +70,18 @@ context("QtlDataset internal helpers") ) } +# `...` reaches QtlDataset(): mafCutoff / keepIndel / keepVariants and the +# other filters are constructor arguments, so a test configures them there +# instead of assigning into the built object's slots. .qh_makeDataset <- function( contexts = c("brain", "liver"), n_samples = 12L, - geno_cov = NULL + geno_cov = NULL, + a1 = rep("A", 6L), + a2 = rep("G", 6L), + ... ) { - gh <- .qh_makeHandle(n_samples = n_samples) + gh <- .qh_makeHandle(n_samples = n_samples, a1 = a1, a2 = a2) pheno <- setNames( lapply(contexts, function(.) .qh_makeSe(n_samples = n_samples)), contexts @@ -80,7 +93,8 @@ context("QtlDataset internal helpers") study = "study1", genotypes = gh, phenotypes = pheno, - genotypeCovariates = geno_cov + genotypeCovariates = geno_cov, + ... ) } @@ -293,7 +307,7 @@ test_that(".qtlResolveVariantRegion: region path expands by cisWindow", { test_that(".qtlVariantIndices: NULL region returns all SNP indices", { qd <- .qh_makeDataset() idx <- pecotmr:::.qtlVariantIndices(qd) - expect_equal(idx, seq_len(nrow(qd@genotypes@snpInfo))) + expect_equal(idx, seq_len(nrow(getSnpInfo(getGenotypeHandle(qd))))) }) test_that(".qtlVariantIndices: filters by chromosome and BP range", { @@ -662,29 +676,29 @@ test_that(".qtlBuildResidualizationDesign: intersects sample sets across blocks" rbinom(n_samples * n_snp, 2, 0.3), nrow = n_samples, ncol = n_snp, - dimnames = list(handle@sampleIds, handle@snpInfo$SNP) + dimnames = list(getSampleIds(handle), getSnpInfo(handle)$SNP) ) sub <- panel[, snpIdx, drop = FALSE] # Build the SE in variants x samples orientation (matches the real impl). rr <- GenomicRanges::GRanges( - seqnames = paste0("chr", handle@snpInfo$CHR[snpIdx]), + seqnames = paste0("chr", getSnpInfo(handle)$CHR[snpIdx]), ranges = IRanges::IRanges( - start = handle@snpInfo$BP[snpIdx], + start = getSnpInfo(handle)$BP[snpIdx], width = 1L ) ) S4Vectors::mcols(rr) <- S4Vectors::DataFrame( - SNP = handle@snpInfo$SNP[snpIdx], - A1 = handle@snpInfo$A1[snpIdx], - A2 = handle@snpInfo$A2[snpIdx] + SNP = getSnpInfo(handle)$SNP[snpIdx], + A1 = getSnpInfo(handle)$A1[snpIdx], + A2 = getSnpInfo(handle)$A2[snpIdx] ) cd <- S4Vectors::DataFrame( - sampleId = handle@sampleIds, - row.names = handle@sampleIds + sampleId = getSampleIds(handle), + row.names = getSampleIds(handle) ) dosage <- t(sub) - rownames(dosage) <- handle@snpInfo$SNP[snpIdx] - colnames(dosage) <- handle@sampleIds + rownames(dosage) <- getSnpInfo(handle)$SNP[snpIdx] + colnames(dosage) <- getSampleIds(handle) SummarizedExperiment::SummarizedExperiment( assays = list(dosage = dosage), rowRanges = rr, @@ -716,8 +730,7 @@ test_that(".qtlExtractBlock: empty snpIdx returns a zero-column block", { }) test_that(".qtlExtractBlock: keepVariants restriction narrows the returned set", { - qd <- .qh_makeDataset() - qd@keepVariants <- c("rs2", "rs4") + qd <- .qh_makeDataset(keepVariants = c("rs2", "rs4")) local_mocked_bindings( extractBlockGenotypes = .qh_mockExtractor(), .package = "pecotmr" @@ -755,8 +768,7 @@ test_that(".qtlExtractBlock: per-call samples arg further narrows the sample set }) test_that(".qtlExtractBlock: keepVariants with empty intersection returns empty block", { - qd <- .qh_makeDataset() - qd@keepVariants <- c("rsGHOST") + qd <- .qh_makeDataset(keepVariants = c("rsGHOST")) local_mocked_bindings( extractBlockGenotypes = .qh_mockExtractor(), .package = "pecotmr" @@ -767,11 +779,12 @@ test_that(".qtlExtractBlock: keepVariants with empty intersection returns empty }) test_that(".qtlExtractBlock: keepIndel = FALSE drops multi-base (indel) variants", { - qd <- .qh_makeDataset() - # Make rs2 (insertion) and rs4 (deletion) indels; the rest stay SNPs. - qd@genotypes@snpInfo$A1[2] <- "AT" - qd@genotypes@snpInfo$A2[4] <- "GC" - qd@keepIndel <- FALSE + # rs2 an insertion, rs4 a deletion; the rest stay SNPs. + qd <- .qh_makeDataset( + a1 = replace(rep("A", 6L), 2L, "AT"), + a2 = replace(rep("G", 6L), 4L, "GC"), + keepIndel = FALSE + ) local_mocked_bindings( extractBlockGenotypes = .qh_mockExtractor(), .package = "pecotmr" @@ -781,8 +794,7 @@ test_that(".qtlExtractBlock: keepIndel = FALSE drops multi-base (indel) variants }) test_that(".qtlExtractBlock: keepIndel = TRUE (default) keeps indel variants", { - qd <- .qh_makeDataset() - qd@genotypes@snpInfo$A1[2] <- "AT" + qd <- .qh_makeDataset(a1 = replace(rep("A", 6L), 2L, "AT")) expect_true(qd@keepIndel) # default local_mocked_bindings( extractBlockGenotypes = .qh_mockExtractor(), @@ -801,11 +813,10 @@ test_that("QtlDataset: keepIndel defaults to TRUE; validity rejects non-scalar", }) test_that(".qtlExtractBlock: mafCutoff drops low-MAF variants", { - qd <- .qh_makeDataset() # The mocked extractor returns binomial(0.3) dosages: realized MAFs hover # around 0.3-0.5 (small sample noise). Cutoff above the realized maximum # drops everything. - qd@mafCutoff <- 0.51 + qd <- .qh_makeDataset(mafCutoff = 0.51) local_mocked_bindings( extractBlockGenotypes = .qh_mockExtractor(), .package = "pecotmr" @@ -815,8 +826,8 @@ test_that(".qtlExtractBlock: mafCutoff drops low-MAF variants", { }) test_that(".qtlExtractBlock: mafCutoff retains variants above the threshold", { - qd <- .qh_makeDataset() - qd@mafCutoff <- 0.4 # realized MAFs include 0.458 and 0.5 + # realized MAFs include 0.458 and 0.5 + qd <- .qh_makeDataset(mafCutoff = 0.4) local_mocked_bindings( extractBlockGenotypes = .qh_mockExtractor(), .package = "pecotmr" @@ -832,35 +843,35 @@ test_that(".qtlExtractBlock: mafCutoff retains variants above the threshold", { function(handle, snpIdx, meanImpute = TRUE) { panel <- matrix( 0, - nrow = length(handle@sampleIds), - ncol = nrow(handle@snpInfo), - dimnames = list(handle@sampleIds, handle@snpInfo$SNP) + nrow = length(getSampleIds(handle)), + ncol = nrow(getSnpInfo(handle)), + dimnames = list(getSampleIds(handle), getSnpInfo(handle)$SNP) ) panel[, "rs1"] <- c(2, 2, 2, 2, 1, 1, 1, 1, 1, 1) # sum 14 -> p = 0.70 panel[, "rs2"] <- c(1, 1, 1, 1, 0, 0, 0, 0, 0, 0) # sum 4 -> p = 0.20 panel[, "rs3"] <- 1 # sum 10 -> p = 0.50 sub <- panel[, snpIdx, drop = FALSE] rr <- GenomicRanges::GRanges( - seqnames = paste0("chr", handle@snpInfo$CHR[snpIdx]), + seqnames = paste0("chr", getSnpInfo(handle)$CHR[snpIdx]), ranges = IRanges::IRanges( - start = handle@snpInfo$BP[snpIdx], + start = getSnpInfo(handle)$BP[snpIdx], width = 1L ) ) S4Vectors::mcols(rr) <- S4Vectors::DataFrame( - SNP = handle@snpInfo$SNP[snpIdx], - A1 = handle@snpInfo$A1[snpIdx], - A2 = handle@snpInfo$A2[snpIdx] + SNP = getSnpInfo(handle)$SNP[snpIdx], + A1 = getSnpInfo(handle)$A1[snpIdx], + A2 = getSnpInfo(handle)$A2[snpIdx] ) dosage <- t(sub) - rownames(dosage) <- handle@snpInfo$SNP[snpIdx] - colnames(dosage) <- handle@sampleIds + rownames(dosage) <- getSnpInfo(handle)$SNP[snpIdx] + colnames(dosage) <- getSampleIds(handle) SummarizedExperiment::SummarizedExperiment( assays = list(dosage = dosage), rowRanges = rr, colData = S4Vectors::DataFrame( - sampleId = handle@sampleIds, - row.names = handle@sampleIds + sampleId = getSampleIds(handle), + row.names = getSampleIds(handle) ) ) } @@ -997,28 +1008,28 @@ context("QtlDataset residualization methods") rbinom(n_samples * n_snp, 2, 0.3), nrow = n_samples, ncol = n_snp, - dimnames = list(handle@sampleIds, handle@snpInfo$SNP) + dimnames = list(getSampleIds(handle), getSnpInfo(handle)$SNP) ) sub <- panel[, snpIdx, drop = FALSE] rr <- GenomicRanges::GRanges( - seqnames = paste0("chr", handle@snpInfo$CHR[snpIdx]), + seqnames = paste0("chr", getSnpInfo(handle)$CHR[snpIdx]), ranges = IRanges::IRanges( - start = handle@snpInfo$BP[snpIdx], + start = getSnpInfo(handle)$BP[snpIdx], width = 1L ) ) S4Vectors::mcols(rr) <- S4Vectors::DataFrame( - SNP = handle@snpInfo$SNP[snpIdx], - A1 = handle@snpInfo$A1[snpIdx], - A2 = handle@snpInfo$A2[snpIdx] + SNP = getSnpInfo(handle)$SNP[snpIdx], + A1 = getSnpInfo(handle)$A1[snpIdx], + A2 = getSnpInfo(handle)$A2[snpIdx] ) cd <- S4Vectors::DataFrame( - sampleId = handle@sampleIds, - row.names = handle@sampleIds + sampleId = getSampleIds(handle), + row.names = getSampleIds(handle) ) dosage <- t(sub) - rownames(dosage) <- handle@snpInfo$SNP[snpIdx] - colnames(dosage) <- handle@sampleIds + rownames(dosage) <- getSnpInfo(handle)$SNP[snpIdx] + colnames(dosage) <- getSampleIds(handle) SummarizedExperiment::SummarizedExperiment( assays = list(dosage = dosage), rowRanges = rr, @@ -2225,9 +2236,9 @@ test_that("multi-range region unions disjoint sub-ranges on one chromosome", { qd <- QtlDataset( study = "S", genotypes = h, - phenotypes = list(ctx = .mr_makeSE(h@sampleIds)) + phenotypes = list(ctx = .mr_makeSE(getSampleIds(h))) ) - bp <- h@snpInfo$BP + bp <- getSnpInfo(h)$BP lo <- min(bp) hi <- max(bp) mid <- lo + (hi - lo) %/% 2L @@ -2257,11 +2268,11 @@ test_that("multi-range region spans chromosomes on a sharded handle", { qd <- QtlDataset( study = "S", genotypes = hs, - phenotypes = list(ctx = .mr_makeSE(hs@sampleIds)) + phenotypes = list(ctx = .mr_makeSE(getSampleIds(hs))) ) - bp <- GenotypeHandle( + bp <- getSnpInfo(GenotypeHandle( plink1Prefix = file.path(test_data_dir, "test_variants") - )@snpInfo$BP + ))$BP lo <- min(bp) hi <- max(bp) r21 <- GenomicRanges::GRanges("chr21", IRanges::IRanges(lo, hi)) @@ -2291,11 +2302,11 @@ test_that("multi-region: single-range extraction is unchanged (regression)", { qd <- QtlDataset( study = "S", genotypes = h, - phenotypes = list(ctx = .mr_makeSE(h@sampleIds)) + phenotypes = list(ctx = .mr_makeSE(getSampleIds(h))) ) - bp <- h@snpInfo$BP + bp <- getSnpInfo(h)$BP r <- GenomicRanges::GRanges("chr21", IRanges::IRanges(min(bp), max(bp))) - expect_equal(.mr_ncol(qd, r), nrow(h@snpInfo)) + expect_equal(.mr_ncol(qd, r), nrow(getSnpInfo(h))) }) test_that(".qtlResolveVariantRegion rejects a non-GRanges / empty region", { @@ -2306,7 +2317,7 @@ test_that(".qtlResolveVariantRegion rejects a non-GRanges / empty region", { qd <- QtlDataset( study = "S", genotypes = h, - phenotypes = list(ctx = .mr_makeSE(h@sampleIds)) + phenotypes = list(ctx = .mr_makeSE(getSampleIds(h))) ) expect_error(getGenotypes(qd, region = "chr21:1-2"), "must be a GRanges") expect_error( @@ -2392,10 +2403,8 @@ test_that(".qtlResolveVariantRegion: region path rejects a non-scalar/negative c # =========================================================================== test_that(".qtlExtractBlock: keepIndel = FALSE with an all-indel panel returns an empty block", { - qd <- .qh_makeDataset() - # Make every variant a multi-base allele so the indel filter drops them all. - qd@genotypes@snpInfo$A1 <- rep("AT", nrow(qd@genotypes@snpInfo)) - qd@keepIndel <- FALSE + # Every variant a multi-base allele, so the indel filter drops them all. + qd <- .qh_makeDataset(a1 = rep("AT", 6L), keepIndel = FALSE) local_mocked_bindings( extractBlockGenotypes = .qh_mockExtractor(), .package = "pecotmr" @@ -2431,47 +2440,47 @@ test_that(".qtlExtractBlock: keepSamples disjoint from the panel returns a zero- # survives the filter (exercising the mean-impute loop). .qh_naExtractor <- function() { function(handle, snpIdx, meanImpute = TRUE) { - ns <- length(handle@sampleIds) - nv <- nrow(handle@snpInfo) + ns <- length(getSampleIds(handle)) + nv <- nrow(getSnpInfo(handle)) set.seed(123L) panel <- matrix( rbinom(ns * nv, 2, 0.4), nrow = ns, ncol = nv, - dimnames = list(handle@sampleIds, handle@snpInfo$SNP) + dimnames = list(getSampleIds(handle), getSnpInfo(handle)$SNP) ) panel["s1", ] <- NA_real_ # fully missing sample -> dropped by imiss panel["s2", "rs2"] <- NA_real_ # scattered NA -> kept then mean-imputed sub <- panel[, snpIdx, drop = FALSE] rr <- GenomicRanges::GRanges( - seqnames = paste0("chr", handle@snpInfo$CHR[snpIdx]), + seqnames = paste0("chr", getSnpInfo(handle)$CHR[snpIdx]), ranges = IRanges::IRanges( - start = handle@snpInfo$BP[snpIdx], + start = getSnpInfo(handle)$BP[snpIdx], width = 1L ) ) S4Vectors::mcols(rr) <- S4Vectors::DataFrame( - SNP = handle@snpInfo$SNP[snpIdx], - A1 = handle@snpInfo$A1[snpIdx], - A2 = handle@snpInfo$A2[snpIdx] + SNP = getSnpInfo(handle)$SNP[snpIdx], + A1 = getSnpInfo(handle)$A1[snpIdx], + A2 = getSnpInfo(handle)$A2[snpIdx] ) dosage <- t(sub) - rownames(dosage) <- handle@snpInfo$SNP[snpIdx] - colnames(dosage) <- handle@sampleIds + rownames(dosage) <- getSnpInfo(handle)$SNP[snpIdx] + colnames(dosage) <- getSampleIds(handle) SummarizedExperiment::SummarizedExperiment( assays = list(dosage = dosage), rowRanges = rr, colData = S4Vectors::DataFrame( - sampleId = handle@sampleIds, - row.names = handle@sampleIds + sampleId = getSampleIds(handle), + row.names = getSampleIds(handle) ) ) } } test_that(".qtlExtractBlock: imissCutoff drops high-missingness samples and mean-imputes the rest", { - qd <- .qh_makeDataset() - qd@imissCutoff <- 0.5 # s1 (100% NA) dropped; s2 (1 NA) kept then imputed + # s1 (100% NA) dropped; s2 (1 NA) kept then imputed + qd <- .qh_makeDataset(imissCutoff = 0.5) local_mocked_bindings( extractBlockGenotypes = .qh_naExtractor(), .package = "pecotmr" @@ -2486,46 +2495,45 @@ test_that(".qtlExtractBlock: imissCutoff drops high-missingness samples and mean # Extractor with one constant (zero-variance) column to drive the xvar filter. .qh_lowVarExtractor <- function() { function(handle, snpIdx, meanImpute = TRUE) { - ns <- length(handle@sampleIds) - nv <- nrow(handle@snpInfo) + ns <- length(getSampleIds(handle)) + nv <- nrow(getSnpInfo(handle)) set.seed(99L) panel <- matrix( rbinom(ns * nv, 2, 0.4), nrow = ns, ncol = nv, - dimnames = list(handle@sampleIds, handle@snpInfo$SNP) + dimnames = list(getSampleIds(handle), getSnpInfo(handle)$SNP) ) panel[, "rs1"] <- 1L # constant column -> variance 0 sub <- panel[, snpIdx, drop = FALSE] rr <- GenomicRanges::GRanges( - seqnames = paste0("chr", handle@snpInfo$CHR[snpIdx]), + seqnames = paste0("chr", getSnpInfo(handle)$CHR[snpIdx]), ranges = IRanges::IRanges( - start = handle@snpInfo$BP[snpIdx], + start = getSnpInfo(handle)$BP[snpIdx], width = 1L ) ) S4Vectors::mcols(rr) <- S4Vectors::DataFrame( - SNP = handle@snpInfo$SNP[snpIdx], - A1 = handle@snpInfo$A1[snpIdx], - A2 = handle@snpInfo$A2[snpIdx] + SNP = getSnpInfo(handle)$SNP[snpIdx], + A1 = getSnpInfo(handle)$A1[snpIdx], + A2 = getSnpInfo(handle)$A2[snpIdx] ) dosage <- t(sub) - rownames(dosage) <- handle@snpInfo$SNP[snpIdx] - colnames(dosage) <- handle@sampleIds + rownames(dosage) <- getSnpInfo(handle)$SNP[snpIdx] + colnames(dosage) <- getSampleIds(handle) SummarizedExperiment::SummarizedExperiment( assays = list(dosage = dosage), rowRanges = rr, colData = S4Vectors::DataFrame( - sampleId = handle@sampleIds, - row.names = handle@sampleIds + sampleId = getSampleIds(handle), + row.names = getSampleIds(handle) ) ) } } test_that(".qtlExtractBlock: xvarCutoff drops near-constant (low-variance) variants", { - qd <- .qh_makeDataset() - qd@xvarCutoff <- 0.01 + qd <- .qh_makeDataset(xvarCutoff = 0.01) local_mocked_bindings( extractBlockGenotypes = .qh_lowVarExtractor(), .package = "pecotmr" @@ -2900,9 +2908,12 @@ test_that("genotype covariates live on the genotype experiment's colData", { cd <- SummarizedExperiment::colData(pecotmr:::.qtlGenotypeSe(qd)) expect_equal(colnames(cd), c("pc1", "pc2")) expect_equal(rownames(cd), getSampleIds(getGenotypeHandle(qd))) - expect_equal(getGenotypeCovariates(qd), pecotmr:::.qtlColDataMatrix( - pecotmr:::.qtlGenotypeSe(qd) - )) + expect_equal( + getGenotypeCovariates(qd), + pecotmr:::.qtlColDataMatrix( + pecotmr:::.qtlGenotypeSe(qd) + ) + ) }) test_that("subsetting by sample keeps the class and its own slots", { @@ -2921,7 +2932,7 @@ test_that("subsetting by context keeps the genotype experiment", { qd <- .qm_makeDataset() # MultiAssayExperiment announces the experiment it drops; here that is # `liver`, which is what was asked for. - expect_warning(out <- qd[, , "brain"], "dropped") + expect_warning(out <- qd[,, "brain"], "dropped") expect_equal( names(MultiAssayExperiment::experiments(out)), c("genotype", "brain") @@ -2933,7 +2944,7 @@ test_that("subsetting by context keeps the genotype experiment", { test_that("subsetting to no context at all is an error", { qd <- .qm_makeDataset() - expect_error(qd[, , "genotype"], "selects no QTL context") + expect_error(qd[,, "genotype"], "selects no QTL context") }) test_that("keepSamples narrows the dataset rather than recording a filter", { @@ -2966,11 +2977,13 @@ test_that("replacing the handle moves the assay's seed with it", { handle <- getGenotypeHandle(qd) handle@path <- "pecotmr://extdata/elsewhere" out <- pecotmr:::.qtlWithGenotypeHandle(qd, handle) - seedPath <- DelayedArray::seed(SummarizedExperiment::assay( - pecotmr:::.qtlGenotypeSe(out), - "dosage" - ))@handle@path - expect_equal(getGenotypeHandle(out)@path, "pecotmr://extdata/elsewhere") + seedPath <- getPath(pecotmr:::.ghSeedHandle( + DelayedArray::seed(SummarizedExperiment::assay( + pecotmr:::.qtlGenotypeSe(out), + "dosage" + )) + )) + expect_equal(getPath(getGenotypeHandle(out)), "pecotmr://extdata/elsewhere") expect_equal(seedPath, "pecotmr://extdata/elsewhere") }) diff --git a/tests/testthat/test_QtlFineMappingResult.R b/tests/testthat/test_QtlFineMappingResult.R index 805af9b3..58351a8c 100644 --- a/tests/testthat/test_QtlFineMappingResult.R +++ b/tests/testthat/test_QtlFineMappingResult.R @@ -14,7 +14,7 @@ test_that("QtlFineMappingResult: builds a collection keyed by 4-tuple", { ) expect_s4_class(res, "QtlFineMappingResult") expect_equal(nrow(res), 2L) - expect_null(res@ldSketch) + expect_null(getLdSketch(res)) }) diff --git a/tests/testthat/test_causalInferencePipeline.R b/tests/testthat/test_causalInferencePipeline.R index 54a7cabe..f236acd6 100644 --- a/tests/testthat/test_causalInferencePipeline.R +++ b/tests/testthat/test_causalInferencePipeline.R @@ -1,7 +1,7 @@ context("causalInferencePipeline") # =========================================================================== -# Strategy: mock extractBlockGenotypes so .cipLdFromSketch returns a real +# Strategy: mock extractBlockGenotypes so .ldFromSketch returns a real # LD matrix on a small panel. Everything else (twasZ, MR, p-value combine) # runs for real on the tiny fixture. # =========================================================================== @@ -29,31 +29,31 @@ context("causalInferencePipeline") function(handle, snpIdx, meanImpute = TRUE) { set.seed(seed) panel <- matrix( - rbinom(n_samples * nrow(handle@snpInfo), 2, 0.3), + rbinom(n_samples * nrow(getSnpInfo(handle)), 2, 0.3), nrow = n_samples, - ncol = nrow(handle@snpInfo), - dimnames = list(handle@sampleIds, handle@snpInfo$SNP) + ncol = nrow(getSnpInfo(handle)), + dimnames = list(getSampleIds(handle), getSnpInfo(handle)$SNP) ) sub <- panel[, snpIdx, drop = FALSE] rr <- GenomicRanges::GRanges( - seqnames = paste0("chr", handle@snpInfo$CHR[snpIdx]), + seqnames = paste0("chr", getSnpInfo(handle)$CHR[snpIdx]), ranges = IRanges::IRanges( - start = handle@snpInfo$BP[snpIdx], + start = getSnpInfo(handle)$BP[snpIdx], width = 1L ) ) S4Vectors::mcols(rr) <- S4Vectors::DataFrame( - SNP = handle@snpInfo$SNP[snpIdx], - A1 = handle@snpInfo$A1[snpIdx], - A2 = handle@snpInfo$A2[snpIdx] + SNP = getSnpInfo(handle)$SNP[snpIdx], + A1 = getSnpInfo(handle)$A1[snpIdx], + A2 = getSnpInfo(handle)$A2[snpIdx] ) cd <- S4Vectors::DataFrame( - sampleId = handle@sampleIds, - row.names = handle@sampleIds + sampleId = getSampleIds(handle), + row.names = getSampleIds(handle) ) dosage <- t(sub) - rownames(dosage) <- handle@snpInfo$SNP[snpIdx] - colnames(dosage) <- handle@sampleIds + rownames(dosage) <- getSnpInfo(handle)$SNP[snpIdx] + colnames(dosage) <- getSampleIds(handle) SummarizedExperiment::SummarizedExperiment( assays = list(dosage = dosage), rowRanges = rr, diff --git a/tests/testthat/test_colocboostPipeline.R b/tests/testthat/test_colocboostPipeline.R index bbb4b159..e154c5e9 100644 --- a/tests/testthat/test_colocboostPipeline.R +++ b/tests/testthat/test_colocboostPipeline.R @@ -112,31 +112,31 @@ context("colocboostPipeline (S4 dispatch)") function(handle, snpIdx, meanImpute = TRUE) { set.seed(seed) panel <- matrix( - rbinom(n_samples * nrow(handle@snpInfo), 2, 0.3), + rbinom(n_samples * nrow(getSnpInfo(handle)), 2, 0.3), nrow = n_samples, - ncol = nrow(handle@snpInfo), - dimnames = list(handle@sampleIds, handle@snpInfo$SNP) + ncol = nrow(getSnpInfo(handle)), + dimnames = list(getSampleIds(handle), getSnpInfo(handle)$SNP) ) sub <- panel[, snpIdx, drop = FALSE] rr <- GenomicRanges::GRanges( - seqnames = paste0("chr", handle@snpInfo$CHR[snpIdx]), + seqnames = paste0("chr", getSnpInfo(handle)$CHR[snpIdx]), ranges = IRanges::IRanges( - start = handle@snpInfo$BP[snpIdx], + start = getSnpInfo(handle)$BP[snpIdx], width = 1L ) ) S4Vectors::mcols(rr) <- S4Vectors::DataFrame( - SNP = handle@snpInfo$SNP[snpIdx], - A1 = handle@snpInfo$A1[snpIdx], - A2 = handle@snpInfo$A2[snpIdx] + SNP = getSnpInfo(handle)$SNP[snpIdx], + A1 = getSnpInfo(handle)$A1[snpIdx], + A2 = getSnpInfo(handle)$A2[snpIdx] ) cd <- S4Vectors::DataFrame( - sampleId = handle@sampleIds, - row.names = handle@sampleIds + sampleId = getSampleIds(handle), + row.names = getSampleIds(handle) ) dosage <- t(sub) - rownames(dosage) <- handle@snpInfo$SNP[snpIdx] - colnames(dosage) <- handle@sampleIds + rownames(dosage) <- getSnpInfo(handle)$SNP[snpIdx] + colnames(dosage) <- getSampleIds(handle) SummarizedExperiment::SummarizedExperiment( assays = list(dosage = dosage), rowRanges = rr, @@ -943,15 +943,15 @@ test_that("colocboostPipeline(MultiStudyQtlDataset): a study with no usable bund # @noRd .cbf_qcd <- function() { data(qtlSumStatsExample, envir = environment()) - suppressWarnings(suppressMessages(summaryStatsQc(qtlSumStatsExample))) + suppressMessages(summaryStatsQc(qtlSumStatsExample)) } # @noRd .cbf_n <- function(ss, ...) { - b <- suppressWarnings(suppressMessages(.cbQtlSumStatsBundle( + b <- suppressMessages(.cbQtlSumStatsBundle( ss, cutoffs = .panelCutoffs(list(...)) - ))) + )) if (length(b) == 0L) 0L else length(b[[1L]]$variantIds) } diff --git a/tests/testthat/test_ctwasPipeline.R b/tests/testthat/test_ctwasPipeline.R index a931c14d..a7fc24be 100644 --- a/tests/testthat/test_ctwasPipeline.R +++ b/tests/testthat/test_ctwasPipeline.R @@ -48,31 +48,31 @@ context("ctwasPipeline") function(handle, snpIdx, meanImpute = TRUE) { set.seed(seed) panel <- matrix( - rbinom(n_samples * nrow(handle@snpInfo), 2, 0.3), + rbinom(n_samples * nrow(getSnpInfo(handle)), 2, 0.3), nrow = n_samples, - ncol = nrow(handle@snpInfo), - dimnames = list(handle@sampleIds, handle@snpInfo$SNP) + ncol = nrow(getSnpInfo(handle)), + dimnames = list(getSampleIds(handle), getSnpInfo(handle)$SNP) ) sub <- panel[, snpIdx, drop = FALSE] rr <- GenomicRanges::GRanges( - seqnames = paste0("chr", handle@snpInfo$CHR[snpIdx]), + seqnames = paste0("chr", getSnpInfo(handle)$CHR[snpIdx]), ranges = IRanges::IRanges( - start = handle@snpInfo$BP[snpIdx], + start = getSnpInfo(handle)$BP[snpIdx], width = 1L ) ) S4Vectors::mcols(rr) <- S4Vectors::DataFrame( - SNP = handle@snpInfo$SNP[snpIdx], - A1 = handle@snpInfo$A1[snpIdx], - A2 = handle@snpInfo$A2[snpIdx] + SNP = getSnpInfo(handle)$SNP[snpIdx], + A1 = getSnpInfo(handle)$A1[snpIdx], + A2 = getSnpInfo(handle)$A2[snpIdx] ) cd <- S4Vectors::DataFrame( - sampleId = handle@sampleIds, - row.names = handle@sampleIds + sampleId = getSampleIds(handle), + row.names = getSampleIds(handle) ) dosage <- t(sub) - rownames(dosage) <- handle@snpInfo$SNP[snpIdx] - colnames(dosage) <- handle@sampleIds + rownames(dosage) <- getSnpInfo(handle)$SNP[snpIdx] + colnames(dosage) <- getSampleIds(handle) SummarizedExperiment::SummarizedExperiment( assays = list(dosage = dosage), rowRanges = rr, @@ -2136,19 +2136,19 @@ test_that("ctwasPipeline: real-engine end-to-end on the bundled example panel", # the previous per-region list dodged this by putting the SAME GWAS block # and the SAME gene under two region ids. geneOne <- twasWeightsRow( - variantIds = gh@snpInfo$SNP[1:5], + variantIds = getSnpInfo(gh)$SNP[1:5], weights = c(0.1, -0.2, 0.05, 0.0, 0.3) ) geneTwo <- twasWeightsRow( - variantIds = gh@snpInfo$SNP[9:13], + variantIds = getSnpInfo(gh)$SNP[9:13], weights = c(0.2, 0.1, -0.15, 0.05, 0.0) ) # A FLAT weight source is placed into blocks by `traitPos`, so each gene # needs the span it sits in. geneSpan <- function(i) { - bp <- gh@snpInfo$BP[i] + bp <- getSnpInfo(gh)$BP[i] GenomicRanges::GRanges( - str_c("chr", gh@snpInfo$CHR[i][[1L]]), + str_c("chr", getSnpInfo(gh)$CHR[i][[1L]]), IRanges::IRanges(min(bp), max(bp)) ) } @@ -2335,7 +2335,7 @@ test_that(".ctwasBuildWeights: maxNumVariants caps the per-gene weight matrix", data(qtlDatasetExample) qd <- qtlDatasetExample gh <- qd@genotypes - vids <- gh@snpInfo$SNP[1:5] + vids <- getSnpInfo(gh)$SNP[1:5] ent <- twasWeightsRow( variantIds = vids, weights = c(0.1, -0.2, 0.05, 0.3, 0.15) @@ -2360,7 +2360,7 @@ test_that(".ctwasBuildWeights: twasWeightCutoff drops low-magnitude variants", { data(qtlDatasetExample) qd <- qtlDatasetExample gh <- qd@genotypes - vids <- gh@snpInfo$SNP[1:5] + vids <- getSnpInfo(gh)$SNP[1:5] ent <- twasWeightsRow( variantIds = vids, # v1 (0.005) and v3 (0.001) will be dropped at cutoff 0.01 diff --git a/tests/testthat/test_fineMappingPipeline.R b/tests/testthat/test_fineMappingPipeline.R index 8e35c135..01877e98 100644 --- a/tests/testthat/test_fineMappingPipeline.R +++ b/tests/testthat/test_fineMappingPipeline.R @@ -34,31 +34,31 @@ context("fineMappingPipeline") function(handle, snpIdx, meanImpute = TRUE) { set.seed(seed) panel <- matrix( - rbinom(n_samples * nrow(handle@snpInfo), 2, 0.3), + rbinom(n_samples * nrow(getSnpInfo(handle)), 2, 0.3), nrow = n_samples, - ncol = nrow(handle@snpInfo), - dimnames = list(handle@sampleIds, handle@snpInfo$SNP) + ncol = nrow(getSnpInfo(handle)), + dimnames = list(getSampleIds(handle), getSnpInfo(handle)$SNP) ) sub <- panel[, snpIdx, drop = FALSE] rr <- GenomicRanges::GRanges( - seqnames = paste0("chr", handle@snpInfo$CHR[snpIdx]), + seqnames = paste0("chr", getSnpInfo(handle)$CHR[snpIdx]), ranges = IRanges::IRanges( - start = handle@snpInfo$BP[snpIdx], + start = getSnpInfo(handle)$BP[snpIdx], width = 1L ) ) S4Vectors::mcols(rr) <- S4Vectors::DataFrame( - SNP = handle@snpInfo$SNP[snpIdx], - A1 = handle@snpInfo$A1[snpIdx], - A2 = handle@snpInfo$A2[snpIdx] + SNP = getSnpInfo(handle)$SNP[snpIdx], + A1 = getSnpInfo(handle)$A1[snpIdx], + A2 = getSnpInfo(handle)$A2[snpIdx] ) cd <- S4Vectors::DataFrame( - sampleId = handle@sampleIds, - row.names = handle@sampleIds + sampleId = getSampleIds(handle), + row.names = getSampleIds(handle) ) dosage <- t(sub) - rownames(dosage) <- handle@snpInfo$SNP[snpIdx] - colnames(dosage) <- handle@sampleIds + rownames(dosage) <- getSnpInfo(handle)$SNP[snpIdx] + colnames(dosage) <- getSampleIds(handle) SummarizedExperiment::SummarizedExperiment( assays = list(dosage = dosage), rowRanges = rr, @@ -705,25 +705,6 @@ test_that(".fmExtractZn: errors on missing SNP / Z / N columns", { expect_error(pecotmr:::.fmExtractZn(gr, "x"), "no N mcol") }) -# =========================================================================== -# .fmLdFromSketch -# =========================================================================== - -test_that(".fmLdFromSketch: returns named LD matrix; missing variants error", { - h <- .fmp_makeHandle() - local_mocked_bindings( - extractBlockGenotypes = .fmp_mockExtractor(), - .package = "pecotmr" - ) - R <- pecotmr:::.fmLdFromSketch(h, c("chr1:100:A:G", "chr1:300:A:G")) - expect_equal(dim(R), c(2L, 2L)) - expect_equal(rownames(R), c("chr1:100:A:G", "chr1:300:A:G")) - expect_error( - pecotmr:::.fmLdFromSketch(h, c("chr1:100:A:G", "ghost")), - "not present in the LD sketch" - ) -}) - # =========================================================================== # fineMappingPipeline(QtlDataset) # =========================================================================== @@ -4691,7 +4672,7 @@ test_that("fineMappingPipeline(QtlDataset): usePCA skips a PC block screened out # @noRd .rssf_qcd <- function() { data(qtlSumStatsExample, envir = environment()) - suppressWarnings(suppressMessages(summaryStatsQc(qtlSumStatsExample))) + suppressMessages(summaryStatsQc(qtlSumStatsExample)) } # @noRd diff --git a/tests/testthat/test_fineMappingRow.R b/tests/testthat/test_fineMappingRow.R index 015e24d2..173abfdc 100644 --- a/tests/testthat/test_fineMappingRow.R +++ b/tests/testthat/test_fineMappingRow.R @@ -1368,8 +1368,11 @@ test_that("every fit credible set gets a row, with its true membership", { alpha = alpha, V = c(0.1, 0.2), lbf_variable = matrix( - c(1, 1, 1, 0, 0, 0, 2, 0, 0, 0), 2, 5, - byrow = TRUE, dimnames = list(c("L1", "L2"), vn) + c(1, 1, 1, 0, 0, 0, 2, 0, 0, 0), + 2, + 5, + byrow = TRUE, + dimnames = list(c("L1", "L2"), vn) ), sets = list( # L2 = {2} sits inside L1 = {1,2,3}: v2 tags to the smaller L2 in @@ -1422,7 +1425,9 @@ test_that("the summary carries the fit's true (gapped) effect index", { ) class(fit) <- "susie" tl <- data.frame( - variant_id = vn, pip = fit$pip, logBF = c(1, 1, 1, 2), + variant_id = vn, + pip = fit$pip, + logBF = c(1, 1, 1, 2), cs_95 = c("susie_2", "susie_2", "susie_2", "susie_1"), stringsAsFactors = FALSE ) @@ -1439,7 +1444,9 @@ test_that("a fit with no stored sets falls back to the per-variant column", { fit <- list(pip = c(0.5, 0.4, 0.1)) class(fit) <- "susie" tl <- data.frame( - variant_id = vn, pip = fit$pip, logBF = c(1, 1, 0), + variant_id = vn, + pip = fit$pip, + logBF = c(1, 1, 0), cs_95 = c("susie_1", "susie_1", "susie_0"), stringsAsFactors = FALSE ) diff --git a/tests/testthat/test_genotypeHandle.R b/tests/testthat/test_genotypeHandle.R index b81a0359..0b41b6af 100644 --- a/tests/testthat/test_genotypeHandle.R +++ b/tests/testthat/test_genotypeHandle.R @@ -61,8 +61,8 @@ test_that("GenotypeHandle: path = .gds uses readGenotypes (gds backend)", { skip_if_not_installed("SNPRelate") h <- GenotypeHandle(path = gds_path) expect_s4_class(h, "GenotypeHandle") - expect_equal(h@format, "gds") - expect_equal(h@nSamples, 100L) + expect_equal(getFormat(h), "gds") + expect_equal(getNSamples(h), 100L) expect_equal(nrow(getSnpInfo(h)), 349L) }) @@ -70,16 +70,16 @@ test_that("GenotypeHandle: path = .vcf.gz uses readGenotypes (vcf backend)", { skip_if_not_installed("VariantAnnotation") h <- GenotypeHandle(path = vcf_path) expect_s4_class(h, "GenotypeHandle") - expect_equal(h@format, "vcf") - expect_equal(h@nSamples, 100L) + expect_equal(getFormat(h), "vcf") + expect_equal(getNSamples(h), 100L) }) test_that("GenotypeHandle: plink1Prefix builds a plink1 handle", { skip_if_not_installed("snpStats") h <- GenotypeHandle(plink1Prefix = plink_prefix) expect_s4_class(h, "GenotypeHandle") - expect_equal(h@format, "plink1") - expect_equal(h@nSamples, 100L) + expect_equal(getFormat(h), "plink1") + expect_equal(getNSamples(h), 100L) expect_equal(nrow(getSnpInfo(h)), 349L) }) @@ -87,8 +87,8 @@ test_that("GenotypeHandle: plink2Prefix builds a plink2 handle", { skip_if_not_installed("pgenlibr") h <- GenotypeHandle(plink2Prefix = plink_prefix) expect_s4_class(h, "GenotypeHandle") - expect_equal(h@format, "plink2") - expect_equal(h@nSamples, 100L) + expect_equal(getFormat(h), "plink2") + expect_equal(getNSamples(h), 100L) }) # =========================================================================== @@ -103,8 +103,8 @@ test_that("GenotypeHandle: bed/bim/fam triplet with matching stems builds plink1 fam = paste0(plink_prefix, ".fam") ) expect_s4_class(h, "GenotypeHandle") - expect_equal(h@format, "plink1") - expect_equal(h@nSamples, 100L) + expect_equal(getFormat(h), "plink1") + expect_equal(getNSamples(h), 100L) }) test_that(".genotypeHandleFromPlink1Triplet: errors when stems disagree", { @@ -137,7 +137,7 @@ test_that("GenotypeHandle: pgen/pvar/psam triplet with matching stems builds pli psam = paste0(plink_prefix, ".psam") ) expect_s4_class(h, "GenotypeHandle") - expect_equal(h@format, "plink2") + expect_equal(getFormat(h), "plink2") }) test_that(".genotypeHandleFromPlink2Triplet: accepts .pvar.zst by stripping the .zst", { @@ -267,7 +267,7 @@ test_that("GenotypeHandle ldMeta: dispatches to gds reader for .gds path", { on.exit(unlink(f), add = TRUE) h <- GenotypeHandle(ldMeta = f, region = "chr21:17513228-17592874") expect_s4_class(h, "GenotypeHandle") - expect_equal(h@format, "gds") + expect_equal(getFormat(h), "gds") }) test_that("GenotypeHandle ldMeta: dispatches to vcf reader for .vcf.gz path", { @@ -276,7 +276,7 @@ test_that("GenotypeHandle ldMeta: dispatches to vcf reader for .vcf.gz path", { on.exit(unlink(f), add = TRUE) h <- GenotypeHandle(ldMeta = f, region = "chr21:17513228-17592874") expect_s4_class(h, "GenotypeHandle") - expect_equal(h@format, "vcf") + expect_equal(getFormat(h), "vcf") }) test_that("GenotypeHandle ldMeta: dispatches to plink1 reader for .bed path", { @@ -285,7 +285,7 @@ test_that("GenotypeHandle ldMeta: dispatches to plink1 reader for .bed path", { on.exit(unlink(f), add = TRUE) h <- GenotypeHandle(ldMeta = f, region = "chr21:17513228-17592874") expect_s4_class(h, "GenotypeHandle") - expect_equal(h@format, "plink1") + expect_equal(getFormat(h), "plink1") }) test_that("GenotypeHandle ldMeta: dispatches to plink2 reader for .pgen path", { @@ -294,7 +294,7 @@ test_that("GenotypeHandle ldMeta: dispatches to plink2 reader for .pgen path", { on.exit(unlink(f), add = TRUE) h <- GenotypeHandle(ldMeta = f, region = "chr21:17513228-17592874") expect_s4_class(h, "GenotypeHandle") - expect_equal(h@format, "plink2") + expect_equal(getFormat(h), "plink2") }) test_that("GenotypeHandle ldMeta: .cor.xz payload is rejected (out of scope)", { @@ -331,7 +331,7 @@ test_that("GenotypeHandle constructs and validates correctly", { pgenPtr = NULL ) expect_s4_class(obj, "GenotypeHandle") - expect_equal(obj@format, "gds") + expect_equal(getFormat(obj), "gds") expect_true(methods::validObject(obj)) }) @@ -383,8 +383,8 @@ test_that("genoMeta (named vector) builds a sharded handle", { ) ) expect_s4_class(h, "GenotypeHandle") - expect_equal(h@format, "plink1") - expect_equal(sort(names(h@chromPaths)), c("21", "22")) + expect_equal(getFormat(h), "plink1") + expect_equal(sort(names(getChromPaths(h))), c("21", "22")) expect_equal(nrow(getSnpInfo(h)), 2L * 349L) }) @@ -398,7 +398,7 @@ test_that("genoMeta meta-file resolves payloads relative to its own directory", metafile ) h <- GenotypeHandle(genoMeta = metafile) - expect_equal(sort(names(h@chromPaths)), c("21", "22")) + expect_equal(sort(names(getChromPaths(h))), c("21", "22")) expect_equal(nrow(getSnpInfo(h)), 2L * 349L) }) @@ -522,11 +522,11 @@ test_that("genoMeta chroms reads only the requested shard", { ) full <- GenotypeHandle(genoMeta = meta) only21 <- GenotypeHandle(genoMeta = meta, chroms = "21") - expect_equal(names(only21@chromPaths), "21") + expect_equal(names(getChromPaths(only21)), "21") expect_equal(nrow(getSnpInfo(only21)), 349L) # chr21 is the first shard, so the kept rows/fileIdx match the full handle's. expect_equal(getSnpInfo(only21)$SNP, getSnpInfo(full)$SNP[1:349]) - expect_equal(only21@snpInfo$fileIdx, full@snpInfo$fileIdx[1:349]) + expect_equal(getSnpInfo(only21)$fileIdx, getSnpInfo(full)$fileIdx[1:349]) }) test_that("genoMeta chroms canonicalises chromosome labels", { @@ -539,7 +539,7 @@ test_that("genoMeta chroms canonicalises chromosome labels", { ), chroms = "chr21" ) - expect_equal(names(h@chromPaths), "21") + expect_equal(names(getChromPaths(h)), "21") }) test_that("genoMeta chroms falls back to all shards when none match", { @@ -553,7 +553,7 @@ test_that("genoMeta chroms falls back to all shards when none match", { ), chroms = "9" ) - expect_equal(sort(names(h@chromPaths)), c("21", "22")) + expect_equal(sort(names(getChromPaths(h))), c("21", "22")) }) test_that("genoMeta chroms skips a shard whose file does not exist", { @@ -565,7 +565,7 @@ test_that("genoMeta chroms skips a shard whose file does not exist", { "21" = "/no/such/chr21/prefix" ) h <- GenotypeHandle(genoMeta = meta, chroms = "22") - expect_equal(names(h@chromPaths), "22") + expect_equal(names(getChromPaths(h)), "22") expect_error(GenotypeHandle(genoMeta = meta, chroms = c("22", "21"))) }) @@ -685,32 +685,32 @@ test_that(".resolveGenotypeShard honours explicit plink1 format and .bed extensi skip_if_not_installed("snpStats") h1 <- pecotmr:::.resolveGenotypeShard(plink_prefix, format = "plink1") expect_s4_class(h1, "GenotypeHandle") - expect_equal(h1@format, "plink1") + expect_equal(getFormat(h1), "plink1") h2 <- pecotmr:::.resolveGenotypeShard(paste0(plink_prefix, ".bed")) - expect_equal(h2@format, "plink1") + expect_equal(getFormat(h2), "plink1") }) test_that(".resolveGenotypeShard honours explicit plink2 format and .pgen extension", { skip_if_not_installed("pgenlibr") h1 <- pecotmr:::.resolveGenotypeShard(plink_prefix, format = "plink2") expect_s4_class(h1, "GenotypeHandle") - expect_equal(h1@format, "plink2") + expect_equal(getFormat(h1), "plink2") h2 <- pecotmr:::.resolveGenotypeShard(paste0(plink_prefix, ".pgen")) - expect_equal(h2@format, "plink2") + expect_equal(getFormat(h2), "plink2") }) test_that(".resolveGenotypeShard dispatches gds via explicit format and .gds extension", { skip_if_not_installed("SNPRelate") h1 <- pecotmr:::.resolveGenotypeShard(gds_path, format = "gds") - expect_equal(h1@format, "gds") + expect_equal(getFormat(h1), "gds") h2 <- pecotmr:::.resolveGenotypeShard(gds_path) - expect_equal(h2@format, "gds") + expect_equal(getFormat(h2), "gds") }) test_that(".resolveGenotypeShard dispatches vcf by extension", { skip_if_not_installed("VariantAnnotation") h <- pecotmr:::.resolveGenotypeShard(vcf_path) - expect_equal(h@format, "vcf") + expect_equal(getFormat(h), "vcf") }) test_that(".resolveGenotypeShard probes the .pgen sidecar for an extension-less prefix", { @@ -796,14 +796,16 @@ test_that("genotypeDelayedArray values match getGenotypes transposed", { test_that("a NULL index means every row or column, not none", { data(qtlDatasetExample, envir = environment()) gh <- getGenotypeHandle(qtlDatasetExample) - seed <- genotypeDelayedArray(gh)@seed + seed <- DelayedArray::seed(genotypeDelayedArray(gh)) full <- S4Arrays::extract_array(seed, list(NULL, NULL)) expect_equal(dim(full), c(nrow(getSnpInfo(gh)), getNSamples(gh))) }) test_that("the seed answers an empty extraction on either axis", { data(qtlDatasetExample, envir = environment()) - seed <- genotypeDelayedArray(getGenotypeHandle(qtlDatasetExample))@seed + seed <- DelayedArray::seed( + genotypeDelayedArray(getGenotypeHandle(qtlDatasetExample)) + ) expect_equal( dim(S4Arrays::extract_array(seed, list(integer(0), NULL))), c(0L, 165L) diff --git a/tests/testthat/test_genotypeIo.R b/tests/testthat/test_genotypeIo.R index ee4fffa4..da47b9cc 100644 --- a/tests/testthat/test_genotypeIo.R +++ b/tests/testthat/test_genotypeIo.R @@ -915,7 +915,7 @@ test_that("matchVariantsToKeep filters to specified variants", { file.path(td, "test_variants"), format = "plink2" ) - vi <- pecotmr:::.snpInfoToVariantInfo(handle@snpInfo) + vi <- pecotmr:::.snpInfoToVariantInfo(getSnpInfo(handle)) # Write a keep file as tab-delimited with chrom/pos columns keep_file <- tempfile(fileext = ".tsv") @@ -938,7 +938,7 @@ test_that("matchVariantsToKeep returns all FALSE for non-matching variants", { file.path(td, "test_variants"), format = "plink2" ) - vi <- pecotmr:::.snpInfoToVariantInfo(handle@snpInfo) + vi <- pecotmr:::.snpInfoToVariantInfo(getSnpInfo(handle)) keep_file <- tempfile(fileext = ".tsv") on.exit(unlink(keep_file), add = TRUE) @@ -1037,7 +1037,7 @@ test_that("matchVariantsToKeep uses position-only matching when no alleles", { file.path(td, "test_variants"), format = "plink2" ) - vi <- pecotmr:::.snpInfoToVariantInfo(handle@snpInfo) + vi <- pecotmr:::.snpInfoToVariantInfo(getSnpInfo(handle)) keep_file <- tempfile(fileext = ".tsv") on.exit(unlink(keep_file), add = TRUE) @@ -1068,7 +1068,7 @@ test_that("readGenotypes loads plink2 handle with all variants", { format = "plink2" ) expect_s4_class(handle, "GenotypeHandle") - expect_equal(handle@nSamples, 100L) + expect_equal(getNSamples(handle), 100L) expect_equal(nrow(getSnpInfo(handle)), 349L) rse <- extractBlockGenotypes(handle, seq_len(nrow(getSnpInfo(handle)))) expect_s4_class(rse, "SummarizedExperiment") @@ -1120,7 +1120,7 @@ test_that("loadGenotypeRegion filters by keep_variants_path for plink2", { file.path(td, "test_variants"), format = "plink2" ) - vi <- pecotmr:::.snpInfoToVariantInfo(handle@snpInfo) + vi <- pecotmr:::.snpInfoToVariantInfo(getSnpInfo(handle)) keep_file <- tempfile(fileext = ".tsv") on.exit(unlink(keep_file), add = TRUE) @@ -1154,8 +1154,8 @@ test_that("readGenotypes plink2 sample names match psam IIDs", { file.path(td, "test_variants"), format = "plink2" ) - expect_true(all(grepl("^(HG|NA)\\d+", handle@sampleIds))) - expect_equal(length(unique(handle@sampleIds)), 100L) + expect_true(all(grepl("^(HG|NA)\\d+", getSampleIds(handle)))) + expect_equal(length(unique(getSampleIds(handle))), 100L) }) # =========================================================================== @@ -1267,9 +1267,9 @@ test_that("readGenotypes creates plink1 handle", { skip_if_not_installed("snpStats") handle <- readGenotypeHandle(plink_prefix, format = "plink1") expect_s4_class(handle, "GenotypeHandle") - expect_equal(handle@format, "plink1") - expect_equal(handle@nSamples, n_samples) - expect_equal(nrow(handle@snpInfo), n_variants) + expect_equal(getFormat(handle), "plink1") + expect_equal(getNSamples(handle), n_samples) + expect_equal(nrow(getSnpInfo(handle)), n_variants) }) test_that("loadGenotypeRegion loads plink1 via dispatch", { @@ -1291,9 +1291,9 @@ test_that("readGenotypes creates plink2 handle", { skip_if_not_installed("pgenlibr") handle <- readGenotypeHandle(plink_prefix, format = "plink2") expect_s4_class(handle, "GenotypeHandle") - expect_equal(handle@format, "plink2") - expect_equal(handle@nSamples, n_samples) - expect_equal(nrow(handle@snpInfo), n_variants) + expect_equal(getFormat(handle), "plink2") + expect_equal(getNSamples(handle), n_samples) + expect_equal(nrow(getSnpInfo(handle)), n_variants) }) test_that("extractBlockGenotypes works for plink2", { @@ -1343,9 +1343,9 @@ test_that("readGenotypes creates vcf handle", { skip_if_not_installed("VariantAnnotation") handle <- readGenotypeHandle(vcf_path, format = "vcf") expect_s4_class(handle, "GenotypeHandle") - expect_equal(handle@format, "vcf") - expect_equal(handle@nSamples, n_samples) - expect_equal(nrow(handle@snpInfo), n_variants) + expect_equal(getFormat(handle), "vcf") + expect_equal(getNSamples(handle), n_samples) + expect_equal(nrow(getSnpInfo(handle)), n_variants) }) test_that("loadGenotypeRegion loads VCF via dispatch", { @@ -1391,9 +1391,9 @@ test_that("readGenotypes creates gds handle", { skip_if_not_installed("gdsfmt") handle <- readGenotypeHandle(gds_path, format = "gds") expect_s4_class(handle, "GenotypeHandle") - expect_equal(handle@format, "gds") - expect_equal(handle@nSamples, n_samples) - expect_equal(nrow(handle@snpInfo), n_variants) + expect_equal(getFormat(handle), "gds") + expect_equal(getNSamples(handle), n_samples) + expect_equal(nrow(getSnpInfo(handle)), n_variants) }) test_that("loadGenotypeRegion loads GDS via dispatch", { @@ -1467,8 +1467,8 @@ test_that("PLINK1 and PLINK2 readGenotypes return consistent alleles", { h1 <- readGenotypeHandle(plink_prefix, format = "plink1") h2 <- readGenotypeHandle(plink_prefix, format = "plink2") - expect_equal(h1@snpInfo$A1, h2@snpInfo$A1) - expect_equal(h1@snpInfo$A2, h2@snpInfo$A2) + expect_equal(getSnpInfo(h1)$A1, getSnpInfo(h2)$A1) + expect_equal(getSnpInfo(h1)$A2, getSnpInfo(h2)$A2) }) # --- loadGenotypeRegion (dispatch) ----------------------------------------- @@ -1534,7 +1534,7 @@ test_that("extractBlockGenotypes returns SummarizedExperiment", { dosage <- SummarizedExperiment::assay(rse, "dosage") # Bioc convention: variants x samples expect_equal(nrow(dosage), min(5L, n_snps)) - expect_equal(ncol(dosage), handle@nSamples) + expect_equal(ncol(dosage), getNSamples(handle)) # rowRanges should have variant info rr <- SummarizedExperiment::rowRanges(rse) expect_true("A1" %in% names(S4Vectors::mcols(rr))) @@ -1653,8 +1653,8 @@ test_that("extractBlockGenotypes returns SummarizedExperiment", { ref22 <- spec$ref("_chr22") shard <- GenotypeHandle(genoMeta = c("21" = spec$p21, "22" = spec$p22)) expect_s4_class(shard, "GenotypeHandle") - expect_equal(shard@format, ref21@format) - expect_equal(sort(names(shard@chromPaths)), c("21", "22")) + expect_equal(getFormat(shard), getFormat(ref21)) + expect_equal(sort(names(getChromPaths(shard))), c("21", "22")) n21 <- nrow(getSnpInfo(ref21)) n22 <- nrow(getSnpInfo(ref22)) expect_equal(nrow(getSnpInfo(shard)), n21 + n22) @@ -1678,7 +1678,7 @@ test_that("extractBlockGenotypes returns SummarizedExperiment", { shard <- GenotypeHandle( genoMeta = c("21" = spec$p21, "22" = spec$p22) ) - n21 <- nrow(ref21@snpInfo) + n21 <- nrow(getSnpInfo(ref21)) em <- extractBlockGenotypes(shard, c(1L, 2L, n21 + 1L, n21 + 2L)) dm <- unname(as.matrix(SummarizedExperiment::assay(em, "dosage"))) expect_equal(nrow(dm), 4L) @@ -1707,8 +1707,8 @@ test_that("single-shard sharded handle equals the single-file handle", { sh <- GenotypeHandle( genoMeta = c("21" = file.path(test_data_dir, "test_variants")) ) - expect_equal(length(sh@chromPaths), 1L) - expect_equal(nrow(sh@snpInfo), nrow(ref@snpInfo)) + expect_equal(length(getChromPaths(sh)), 1L) + expect_equal(nrow(getSnpInfo(sh)), nrow(getSnpInfo(ref))) expect_equal(.shardDose(sh, 1:10), .shardDose(ref, 1:10)) }) @@ -1731,8 +1731,11 @@ test_that("genoMeta meta-file form matches the named-vector form", { "22" = file.path(td_abs, "test_variants_chr22") ) ) - expect_equal(nrow(hFile@snpInfo), nrow(hVec@snpInfo)) - expect_equal(sort(names(hFile@chromPaths)), sort(names(hVec@chromPaths))) + expect_equal(nrow(getSnpInfo(hFile)), nrow(getSnpInfo(hVec))) + expect_equal( + sort(names(getChromPaths(hFile))), + sort(names(getChromPaths(hVec))) + ) expect_equal(.shardDose(hFile, 1:5), .shardDose(hVec, 1:5)) }) @@ -1815,8 +1818,8 @@ test_that("extractBlockGenotypes on a sharded handle handles an empty block", { expect_s4_class(se, "SummarizedExperiment") dosage <- SummarizedExperiment::assay(se, "dosage") expect_equal(nrow(dosage), 0L) - expect_equal(ncol(dosage), shard@nSamples) - expect_equal(colnames(dosage), shard@sampleIds) + expect_equal(ncol(dosage), getNSamples(shard)) + expect_equal(colnames(dosage), getSampleIds(shard)) }) test_that("sharded extraction errors for a chromosome with no payload", { @@ -1830,8 +1833,8 @@ test_that("sharded extraction errors for a chromosome with no payload", { ) # Drop the chr22 payload but keep its variants in @snpInfo, so routing a # chr22 request finds no per-chromosome file. - shard@chromPaths <- shard@chromPaths["21"] - idx22 <- which(pecotmr:::canonChrom(shard@snpInfo$CHR) == "22")[1] + shard@chromPaths <- getChromPaths(shard)["21"] + idx22 <- which(pecotmr:::canonChrom(getSnpInfo(shard)$CHR) == "22")[1] expect_error( extractBlockGenotypes(shard, idx22), "no per-chromosome file for chromosome" @@ -2196,7 +2199,7 @@ test_that("loadGenotypeRegion warns on non-integer dosages without a sidecar", { # non-integer values and emits its warning. testthat::local_mocked_bindings( .extractBlockPlink2 = function(handle, snpIdx) { - matrix(0.5, nrow = handle@nSamples, ncol = length(snpIdx)) + matrix(0.5, nrow = getNSamples(handle), ncol = length(snpIdx)) }, .package = "pecotmr" ) @@ -2341,8 +2344,12 @@ test_that("onDisk LD matches the in-memory path on complete data", { format = "gds" ) idx <- c(4L, 1L, 6L, 2L, 5L, 3L) - onDisk <- computeLd(handle, snpIdx = idx, backend = "snprelate", - onDisk = TRUE) + onDisk <- computeLd( + handle, + snpIdx = idx, + backend = "snprelate", + onDisk = TRUE + ) inMem <- computeLd(handle, snpIdx = idx) # The fixture has no missing calls, which is the only regime where the # two are meant to agree. @@ -2356,8 +2363,7 @@ test_that("onDisk LD labels and orders by the request", { format = "gds" ) idx <- c(4L, 1L, 6L, 2L, 5L, 3L) - R <- computeLd(handle, snpIdx = idx, backend = "snprelate", - onDisk = TRUE) + R <- computeLd(handle, snpIdx = idx, backend = "snprelate", onDisk = TRUE) # snpgdsLDMat() returns file order and labels nothing; both are corrected. expect_identical(rownames(R), getSnpInfo(handle)$SNP[idx]) expect_identical(colnames(R), rownames(R)) @@ -2434,9 +2440,14 @@ test_that("readGenotypes returns a panel, not a handle", { ) # variants x samples, the Bioconductor orientation. expect_equal(nrow(panel), nrow(getSnpInfo(handle))) - expect_equal(ncol(panel), length(pecotmr:::.ghSeedHandle( - DelayedArray::seed(SummarizedExperiment::assay(panel, "dosage")) - )@sampleIds)) + expect_equal( + ncol(panel), + length(getSampleIds( + pecotmr:::.ghSeedHandle( + DelayedArray::seed(SummarizedExperiment::assay(panel, "dosage")) + ) + )) + ) }) test_that("the panel's assay is delayed and its seed is the handle", { diff --git a/tests/testthat/test_h2EstimationWrappers.R b/tests/testthat/test_h2EstimationWrappers.R index 10ba1450..da310dc5 100644 --- a/tests/testthat/test_h2EstimationWrappers.R +++ b/tests/testthat/test_h2EstimationWrappers.R @@ -468,13 +468,13 @@ test_that("estimateh2 with method='lder' returns H2Estimate with correct slots", result <- estimateH2(ss, eigen_ref, method = "lder") expect_s4_class(result, "H2Estimate") - expect_true(is.numeric(result@h2)) - expect_true(is.numeric(result@h2Se)) - expect_true(is.numeric(result@intercept)) - expect_true(is.numeric(result@interceptSe)) - expect_equal(result@method, "lder") - expect_equal(result@nSnps, nSnps(ss)) - expect_equal(result@traitName, "test") + expect_true(is.numeric(getH2(result))) + expect_true(is.numeric(getH2Se(result))) + expect_true(is.numeric(getIntercept(result))) + expect_true(is.numeric(getInterceptSe(result))) + expect_equal(getMethodNames(result), "lder") + expect_equal(getNSnps(result), nSnps(ss)) + expect_equal(getTraitName(result), "test") }) test_that("estimateh2 with var_y correction runs without error", { @@ -488,7 +488,7 @@ test_that("estimateh2 with var_y correction runs without error", { result <- estimateH2(ss, eigen_ref, method = "lder") expect_s4_class(result, "H2Estimate") - expect_true(is.numeric(result@h2)) + expect_true(is.numeric(getH2(result))) }) test_that("estimateh2 with method='gldsc' returns H2Estimate", { @@ -497,11 +497,11 @@ test_that("estimateh2 with method='gldsc' returns H2Estimate", { result <- estimateH2(ss, score_ref, method = "gldsc") expect_s4_class(result, "H2Estimate") - expect_true(is.numeric(result@h2)) - expect_true(is.numeric(result@h2Se)) - expect_equal(result@method, "gldsc") - expect_equal(result@nSnps, nSnps(ss)) - expect_equal(result@traitName, "test") + expect_true(is.numeric(getH2(result))) + expect_true(is.numeric(getH2Se(result))) + expect_equal(getMethodNames(result), "gldsc") + expect_equal(getNSnps(result), nSnps(ss)) + expect_equal(getTraitName(result), "test") }) test_that("estimateh2 with method='hdl' returns H2Estimate", { @@ -514,11 +514,11 @@ test_that("estimateh2 with method='hdl' returns H2Estimate", { ) expect_s4_class(result, "H2Estimate") - expect_true(is.numeric(result@h2)) - expect_true(is.numeric(result@h2Se)) - expect_equal(result@method, "hdl") - expect_equal(result@nSnps, nSnps(ss)) - expect_equal(result@traitName, "test") + expect_true(is.numeric(getH2(result))) + expect_true(is.numeric(getH2Se(result))) + expect_equal(getMethodNames(result), "hdl") + expect_equal(getNSnps(result), nSnps(ss)) + expect_equal(getTraitName(result), "test") }) # =========================================================================== @@ -1871,7 +1871,7 @@ test_that("gldscUnivariate with baseline-only annotations yields NULL scoreStats test_that("h2EstimateToSldscTrait assigns category names to unnamed tauBlocks", { h2_obj <- make_test_h2estimate(with_enrichment = TRUE) - tb <- h2_obj@tauBlocks + tb <- getTauBlocks(h2_obj) colnames(tb) <- NULL h2_obj@tauBlocks <- tb result <- h2EstimateToSldscTrait(h2_obj) @@ -1890,7 +1890,9 @@ test_that("estimateH2 errors when study is omitted for a multi-study collection" df <- data.frame( SNP = names(eigen_ref), CHR = sub( - "^chr", "", as.character(GenomicRanges::seqnames(eigen_ref)) + "^chr", + "", + as.character(GenomicRanges::seqnames(eigen_ref)) ), BP = GenomicRanges::start(eigen_ref), A1 = as.character(S4Vectors::mcols(eigen_ref)$A1), diff --git a/tests/testthat/test_jointEngine.R b/tests/testthat/test_jointEngine.R index dc9a2158..f120d333 100644 --- a/tests/testthat/test_jointEngine.R +++ b/tests/testthat/test_jointEngine.R @@ -443,11 +443,6 @@ test_that(".runJointCell: cross-context twas CV-only rows (fitFullData=FALSE)", paste0("c", seq_len(k)) ) ) - R <- diag(p) - dimnames(R) <- list( - sprintf("chr1:%d:A:G", 100L * (seq_len(p))), - sprintf("chr1:%d:A:G", 100L * (seq_len(p))) - ) new( "SumStatsJointGroup", conditions = data.frame( @@ -457,11 +452,20 @@ test_that(".runJointCell: cross-context twas CV-only rows (fitFullData=FALSE)", stringsAsFactors = FALSE ), Z = Z, - R = R, + # The group carries the LD REFERENCE now; tests that fit mock + # .ldFromSketch to supply the matrix (see .je_mockSsLd). + ldSketch = NULL, N = c(100, 120) ) } +# Identity LD over whatever variants are asked for -- the shape .je_ssGroup +# used to hard-code into its R slot. +.je_mockSsLd <- function(sketch, vids, ...) { + diag(length(vids)) |> + `dimnames<-`(list(vids, vids)) +} + .je_ssCell <- function(groups) { new( "JointDispatchCell", @@ -473,6 +477,11 @@ test_that(".runJointCell: cross-context twas CV-only rows (fitFullData=FALSE)", } test_that(".runJointCell: cross-context FM sumstats (mvsusie_rss) -> per-context", { + # The group carries the LD reference; supply the matrix it derives. + local_mocked_bindings( + .ldFromSketch = .je_mockSsLd, + .package = "pecotmr" + ) set.seed(6) cell <- .je_ssCell(list(.je_mkSsGroup("G1"))) pipe <- new("FmJointPipeline", config = list(coverage = 0.95)) @@ -505,6 +514,11 @@ test_that(".runJointCell: cross-context FM sumstats (mvsusie_rss) -> per-context }) test_that(".runJointCell: cross-context twas sumstats (mr.mash.rss) -> per-context", { + # The group carries the LD reference; supply the matrix it derives. + local_mocked_bindings( + .ldFromSketch = .je_mockSsLd, + .package = "pecotmr" + ) set.seed(7) cell <- .je_ssCell(list(.je_mkSsGroup("G1"))) pipe <- new("TwasJointPipeline", config = list()) @@ -538,6 +552,11 @@ test_that(".runJointCell: cross-context twas sumstats (mr.mash.rss) -> per-conte }) test_that("fitJointGroup(SumStats, Twas): real mr.mash-rss keys stat$n (regression)", { + # The group carries the LD reference; supply the matrix it derives. + local_mocked_bindings( + .ldFromSketch = .je_mockSsLd, + .package = "pecotmr" + ) skip_if_not_installed("mr.mashr") # The sibling tests above mock mrmashRssWeights, so the real reader # (.mrmashRssStats, which keys on stat$n) was never exercised. The method @@ -565,7 +584,7 @@ test_that("fitJointGroup(SumStats, Twas): real mr.mash-rss keys stat$n (regressi stringsAsFactors = FALSE ), Z = Z, - R = R, + ldSketch = NULL, N = rep(nObs, K) ) pipe <- new( @@ -605,6 +624,11 @@ test_that(".lookupJointCell: present cells resolve, absent cells error", { # ---- cross-study pattern (study jointed; sumstats-only) --------------------- test_that(".runJointCell: cross-study (twas sumstats) -> per-study rows + jointStudies", { + # The group carries the LD reference; supply the matrix it derives. + local_mocked_bindings( + .ldFromSketch = .je_mockSsLd, + .package = "pecotmr" + ) set.seed(10) Z <- matrix( rnorm(6), @@ -626,7 +650,7 @@ test_that(".runJointCell: cross-study (twas sumstats) -> per-study rows + jointS stringsAsFactors = FALSE ), Z = Z, - R = R, + ldSketch = NULL, N = c(100, 120) ) cell <- new( @@ -849,7 +873,7 @@ test_that("fitJointGroup(Individual, Fm): fsusie returns one entry per trait", { ), X = X, Y = Y, - pos = c(100, 200) + traitPos = c(100, 200) ) pipe <- new("FmJointPipeline", config = list(coverage = 0.95)) captured <- NULL @@ -897,6 +921,11 @@ test_that("fitJointGroup(Individual, Fm): fsusie without pos errors; unknown tok }) test_that(".runJointCell: composed/sumstats (context+trait vary) -> per-tuple rows", { + # The group carries the LD reference; supply the matrix it derives. + local_mocked_bindings( + .ldFromSketch = .je_mockSsLd, + .package = "pecotmr" + ) set.seed(13) Z <- matrix( rnorm(9), @@ -922,7 +951,7 @@ test_that(".runJointCell: composed/sumstats (context+trait vary) -> per-tuple ro "SumStatsJointGroup", conditions = conds, Z = Z, - R = R, + ldSketch = NULL, N = c(100, 100, 120) ) cell <- new( @@ -995,8 +1024,8 @@ test_that(".runJointCell: composed/sumstats (context+trait vary) -> per-tuple ro # One method's per-condition entries; CV predictions correlate with Y by predCor # so the R^2 the layer reads is controllable. .je_ensEntries <- function(group, predCor) { - Y <- group@Y - vars <- colnames(group@X) + Y <- .jgY(group) + vars <- colnames(.jgX(group)) lapply(seq_len(ncol(Y)), function(r) { pr <- predCor * Y[, r] + rnorm(nrow(Y), sd = 0.3) names(pr) <- rownames(Y) @@ -1146,7 +1175,7 @@ test_that("fitJointGroup(twas): FM-derived method reuses fine-mapping's CV (hand "TwasJointPipeline", config = list(cvFolds = 2L, ensemble = FALSE) ) - samp <- rownames(g@X) + samp <- rownames(.jgX(g)) fmCv <- list( samplePartition = data.frame( Sample = samp, @@ -1249,7 +1278,7 @@ test_that("fitJointGroup(twas): FM-derived method reuses fine-mapping's CV (hand variantIds = sprintf("chr1:%d:A:G", 100L * (seq_len(p))) ) } -.je_mockLd <- function(sketch, vids) { +.je_mockLd <- function(sketch, vids, ...) { matrix(0, length(vids), length(vids), dimnames = list(vids, vids)) } @@ -1318,8 +1347,8 @@ test_that(".enumCrossContextIndividual: one group per trait in >= 2 contexts", { g <- pecotmr:::.enumCrossContextIndividual(NULL, scope) expect_length(g, 1L) # only G1 survives expect_s4_class(g[[1L]], "IndividualJointGroup") - expect_equal(as.character(g[[1L]]@conditions$context), c("c1", "c2")) - expect_equal(as.character(g[[1L]]@conditions$trait), c("G1", "G1")) + expect_equal(as.character(.jgConditions(g[[1L]])$context), c("c1", "c2")) + expect_equal(as.character(.jgConditions(g[[1L]])$trait), c("G1", "G1")) }) test_that(".enumCrossContextIndividual: study not in scope / < 2 contexts -> empty", { @@ -1352,22 +1381,22 @@ test_that(".enumCrossContextSumstats: groups per (study, trait) with >= 2 contex traits = list(S = "t1") ) local_mocked_bindings( - getLdSketch = function(x) "SKETCH", + getLdSketch = function(x) NULL, .buildJointSumstatZMatrix = .je_mockJointZ, - .fmLdFromSketch = .je_mockLd, + .ldFromSketch = .je_mockLd, .package = "pecotmr" ) g <- pecotmr:::.enumCrossContextSumstats(df, scope) expect_length(g, 1L) expect_s4_class(g[[1L]], "SumStatsJointGroup") - expect_equal(as.character(g[[1L]]@conditions$context), c("c1", "c2")) + expect_equal(as.character(.jgConditions(g[[1L]])$context), c("c1", "c2")) }) test_that(".enumCrossContextSumstats: < 2 contexts and < 2 tuple rows skip", { local_mocked_bindings( - getLdSketch = function(x) "SKETCH", + getLdSketch = function(x) NULL, .buildJointSumstatZMatrix = .je_mockJointZ, - .fmLdFromSketch = .je_mockLd, + .ldFromSketch = .je_mockLd, .package = "pecotmr" ) # study scoped to two contexts but only one row present -> < 2 tupleRows skip @@ -1429,9 +1458,10 @@ test_that(".enumCrossTraitIndividual: one group per context with >= 2 traits + p ) g <- pecotmr:::.enumCrossTraitIndividual(NULL, scope) expect_length(g, 1L) - expect_equal(as.character(g[[1L]]@conditions$context), c("c1", "c1")) - expect_equal(as.character(g[[1L]]@conditions$trait), c("G1", "G2")) - expect_equal(g[[1L]]@pos, c(100, 200)) # rowRanges midpoints + expect_equal(as.character(.jgConditions(g[[1L]])$context), c("c1", "c1")) + expect_equal(as.character(.jgConditions(g[[1L]])$trait), c("G1", "G2")) + # one per TRAIT (rowRanges midpoints), not per variant + expect_equal(.jgTraitPos(g[[1L]]), c(100, 200)) }) test_that(".enumCrossTraitIndividual: study not in scope -> empty", { @@ -1453,14 +1483,14 @@ test_that(".enumCrossTraitSumstats: groups per (study, context) with >= 2 traits traits = list(S = c("t1", "t2")) ) local_mocked_bindings( - getLdSketch = function(x) "SKETCH", + getLdSketch = function(x) NULL, .buildJointSumstatZMatrix = .je_mockJointZ, - .fmLdFromSketch = .je_mockLd, + .ldFromSketch = .je_mockLd, .package = "pecotmr" ) g <- pecotmr:::.enumCrossTraitSumstats(df, scope) expect_length(g, 1L) - expect_equal(as.character(g[[1L]]@conditions$trait), c("t1", "t2")) + expect_equal(as.character(.jgConditions(g[[1L]])$trait), c("t1", "t2")) # < 2 traits present -> skip df1 <- .je_ssDf(studies = "S", contexts = "c1", traits = "t1") expect_length(pecotmr:::.enumCrossTraitSumstats(df1, scope), 0L) @@ -1474,14 +1504,14 @@ test_that(".enumCrossStudySumstats: group per (context, trait) in >= 2 studies", traits = list(S1 = "t1", S2 = "t1") ) local_mocked_bindings( - getLdSketch = function(x) "SKETCH", + getLdSketch = function(x) NULL, .buildJointSumstatZMatrix = .je_mockJointZ, - .fmLdFromSketch = .je_mockLd, + .ldFromSketch = .je_mockLd, .package = "pecotmr" ) g <- pecotmr:::.enumCrossStudySumstats(df, scope) expect_length(g, 1L) - expect_equal(as.character(g[[1L]]@conditions$study), c("S1", "S2")) + expect_equal(as.character(.jgConditions(g[[1L]])$study), c("S1", "S2")) # only one study in scope for the tuple -> filtered to < 2 -> skip scope1 <- list( studies = c("S1", "S2"), @@ -1529,8 +1559,14 @@ test_that(".enumComposedIndividual: one group joining every (context, trait) tup ) g <- pecotmr:::.enumComposedIndividual(NULL, scope) expect_length(g, 1L) - expect_equal(as.character(g[[1L]]@conditions$context), c("c1", "c1", "c2")) - expect_equal(as.character(g[[1L]]@conditions$trait), c("gA", "gB", "gA")) + expect_equal( + as.character(.jgConditions(g[[1L]])$context), + c("c1", "c1", "c2") + ) + expect_equal( + as.character(.jgConditions(g[[1L]])$trait), + c("gA", "gB", "gA") + ) }) test_that(".enumComposedIndividual: study not in scope / NULL xy -> empty", { @@ -1591,7 +1627,7 @@ test_that(".enumUnivariateIndividual: one 1-condition group per (context, trait) g <- pecotmr:::.enumUnivariateIndividual(NULL, scope) expect_length(g, 4L) # 2 ctx x 2 traits expect_true(all( - vapply(g, function(x) nrow(x@conditions), integer(1)) == 1L + vapply(g, function(x) nrow(.jgConditions(x)), integer(1)) == 1L )) }) @@ -1645,7 +1681,7 @@ test_that(".enumComposedSumstats: one group per fixed-axis row block", { traits = list(S = "t1") ) local_mocked_bindings( - getLdSketch = function(x) "SKETCH", + getLdSketch = function(x) NULL, .enumerateComposedSumstatGroups = function(spec, data, scope) { list( groups = list(c(1L, 2L)), @@ -1655,7 +1691,7 @@ test_that(".enumComposedSumstats: one group per fixed-axis row block", { ) }, .buildJointSumstatZMatrix = .je_mockJointZ, - .fmLdFromSketch = .je_mockLd, + .ldFromSketch = .je_mockLd, .package = "pecotmr" ) g <- pecotmr:::.enumComposedSumstats( @@ -1664,12 +1700,12 @@ test_that(".enumComposedSumstats: one group per fixed-axis row block", { args = list(axes = c("context", "trait")) ) expect_length(g, 1L) - expect_equal(as.character(g[[1L]]@conditions$context), c("c1", "c2")) + expect_equal(as.character(.jgConditions(g[[1L]])$context), c("c1", "c2")) }) test_that(".enumComposedSumstats: NULL group index and singleton blocks skip", { local_mocked_bindings( - getLdSketch = function(x) "SKETCH", + getLdSketch = function(x) NULL, .enumerateComposedSumstatGroups = function(spec, data, scope) NULL, .package = "pecotmr" ) @@ -1681,7 +1717,7 @@ test_that(".enumComposedSumstats: NULL group index and singleton blocks skip", { 0L ) # gi NULL (646) local_mocked_bindings( - getLdSketch = function(x) "SKETCH", + getLdSketch = function(x) NULL, .enumerateComposedSumstatGroups = function(spec, data, scope) { list( groups = list(1L), @@ -1691,7 +1727,7 @@ test_that(".enumComposedSumstats: NULL group index and singleton blocks skip", { ) }, .buildJointSumstatZMatrix = .je_mockJointZ, - .fmLdFromSketch = .je_mockLd, + .ldFromSketch = .je_mockLd, .package = "pecotmr" ) expect_length( @@ -1733,7 +1769,7 @@ test_that("fitJointGroup(Individual, Fm): fsusie honest per-fold CV is attached" ), X = X, Y = Y, - pos = c(100, 200) + traitPos = c(100, 200) ) pipe <- new("FmJointPipeline", config = list(coverage = 0.95, cvFolds = 3)) cvCalled <- FALSE @@ -1830,6 +1866,11 @@ test_that("fitJointGroup(SumStats, Fm): fsusie and unknown tokens error", { }) test_that("fitJointGroup(SumStats, Fm): a reweighted-prior residual variance is threaded", { + # The group carries the LD reference; supply the matrix it derives. + local_mocked_bindings( + .ldFromSketch = .je_mockSsLd, + .package = "pecotmr" + ) grp <- .je_mkSsGroup("G1") pipe <- new("FmJointPipeline", config = list(coverage = 0.95)) captured <- NULL @@ -1904,6 +1945,11 @@ test_that("fitJointGroup(twas): spike-and-slab pi feeds bayes_b probIn", { }) test_that("fitJointGroup(SumStats, twas): a vector weight without rownames falls back to Z rows", { + # The group carries the LD reference; supply the matrix it derives. + local_mocked_bindings( + .ldFromSketch = .je_mockSsLd, + .package = "pecotmr" + ) grp <- .je_mkSsGroup("G1", p = 3L, k = 2L) pipe <- new("TwasJointPipeline", config = list()) local_mocked_bindings( @@ -1925,7 +1971,7 @@ test_that("fitJointGroup(SumStats, twas): a vector weight without rownames falls expect_s4_class(res, "TwasWeights") expect_equal( getVariantIds(pecotmr:::.collectionEntry(res, 1L)), - rownames(grp@Z) + rownames(.jgZ(grp)) ) # fallback vids }) @@ -1955,7 +2001,7 @@ test_that(".runJointCell: a fitter returning all-NULL entries yields no rows -> # fitJointGroup returns a list of NULLs (every condition screened out). local_mocked_bindings( fitJointGroup = function(group, pipeline, token, args) { - vector("list", nrow(group@conditions)) + vector("list", nrow(.jgConditions(group))) }, .package = "pecotmr" ) @@ -2000,8 +2046,8 @@ test_that(".twasEnsembleLayer: entries lacking CV predictions are skipped", { good <- .je_ensEntries(g, 0.85) # A method whose entry carries no CV predictions -> contributes nothing. noCv <- list(twasWeightsRow( - variantIds = colnames(g@X), - weights = rnorm(ncol(g@X)), + variantIds = colnames(.jgX(g)), + weights = rnorm(ncol(.jgX(g))), cvResult = NULL )) ens <- pecotmr:::.twasEnsembleLayer( @@ -2052,7 +2098,7 @@ test_that(".twasEnsembleLayer: unnamed ensemble weights fall back to a method's alpha ) { list( - ensembleTwasWeights = as.numeric(rep(0.1, ncol(g@X))), # no names + ensembleTwasWeights = as.numeric(rep(0.1, ncol(.jgX(g)))), # no names methodCoef = c(0.5, 0.5), methodPerformance = c(0.8, 0.7) ) @@ -2070,7 +2116,7 @@ test_that(".twasEnsembleLayer: unnamed ensemble weights fall back to a method's ) ) expect_s4_class(ens[[1L]], "TwasWeightsRow") - expect_equal(getVariantIds(ens[[1L]]), colnames(g@X)) # fallback ids + expect_equal(getVariantIds(ens[[1L]]), colnames(.jgX(g))) # fallback ids }) # ============================================================================= @@ -2156,8 +2202,8 @@ test_that(".twasGroupArgs: takes the CV partition from the fine-mapping CV when g <- .je_mkGroup("G1") # IndividualJointGroup pipe <- new("TwasJointPipeline", config = list(cvFolds = 2L)) sp <- data.frame( - Sample = rownames(g@X), - Fold = rep(1:2, length.out = nrow(g@X)), + Sample = rownames(.jgX(g)), + Fold = rep(1:2, length.out = nrow(.jgX(g))), stringsAsFactors = FALSE ) local_mocked_bindings( diff --git a/tests/testthat/test_jointSpecification.R b/tests/testthat/test_jointSpecification.R index 3a6b8d18..e98ee421 100644 --- a/tests/testthat/test_jointSpecification.R +++ b/tests/testthat/test_jointSpecification.R @@ -735,31 +735,31 @@ context("joint dispatchers (fineMappingDispatcher / twasDispatcher)") function(handle, snpIdx, meanImpute = TRUE) { set.seed(seed) panel <- matrix( - rbinom(n_samples * nrow(handle@snpInfo), 2, 0.3), + rbinom(n_samples * nrow(getSnpInfo(handle)), 2, 0.3), nrow = n_samples, - ncol = nrow(handle@snpInfo), - dimnames = list(handle@sampleIds, handle@snpInfo$SNP) + ncol = nrow(getSnpInfo(handle)), + dimnames = list(getSampleIds(handle), getSnpInfo(handle)$SNP) ) sub <- panel[, snpIdx, drop = FALSE] rr <- GenomicRanges::GRanges( - seqnames = paste0("chr", handle@snpInfo$CHR[snpIdx]), + seqnames = paste0("chr", getSnpInfo(handle)$CHR[snpIdx]), ranges = IRanges::IRanges( - start = handle@snpInfo$BP[snpIdx], + start = getSnpInfo(handle)$BP[snpIdx], width = 1L ) ) S4Vectors::mcols(rr) <- S4Vectors::DataFrame( - SNP = handle@snpInfo$SNP[snpIdx], - A1 = handle@snpInfo$A1[snpIdx], - A2 = handle@snpInfo$A2[snpIdx] + SNP = getSnpInfo(handle)$SNP[snpIdx], + A1 = getSnpInfo(handle)$A1[snpIdx], + A2 = getSnpInfo(handle)$A2[snpIdx] ) cd <- S4Vectors::DataFrame( - sampleId = handle@sampleIds, - row.names = handle@sampleIds + sampleId = getSampleIds(handle), + row.names = getSampleIds(handle) ) dosage <- t(sub) - rownames(dosage) <- handle@snpInfo$SNP[snpIdx] - colnames(dosage) <- handle@sampleIds + rownames(dosage) <- getSnpInfo(handle)$SNP[snpIdx] + colnames(dosage) <- getSampleIds(handle) SummarizedExperiment::SummarizedExperiment( assays = list(dosage = dosage), rowRanges = rr, diff --git a/tests/testthat/test_ld.R b/tests/testthat/test_ld.R index d1451dc1..cf059f71 100644 --- a/tests/testthat/test_ld.R +++ b/tests/testthat/test_ld.R @@ -342,10 +342,10 @@ test_that("partitionLdMatrix validates block structure properly", { # Rebuild LdData with modified matrix and block metadata invalid_ld_data <- new( "LdData", + getVariantInfo(ld_data), correlation = ldmat, genotypeHandle = NULL, - variants = ld_data@variants, - snpIdx = ld_data@snpIdx, + snpIdx = getSnpIdx(ld_data), blockMetadata = bm ) @@ -412,7 +412,7 @@ test_that("partitionLdMatrix handles row/column name mismatches", { colnames(ldmat) <- NULL mismatched_ld_data <- LdData( correlation = ldmat, - variants = ld_data@variants, + variants = getVariantInfo(ld_data), blockMetadata = getBlockMetadata(ld_data) ) @@ -2755,9 +2755,9 @@ test_that("loadLdSketch: returns LdData with raw genotypes and metadata", { # Store genotype matrix directly in genotype_handle (matching loadLdSketch output) mock_ld_data <- new( "LdData", + variants_gr, correlation = NULL, genotypeHandle = X, - variants = variants_gr, snpIdx = NULL, blockMetadata = blockMetadata ) @@ -2821,9 +2821,9 @@ test_that("loadLdSketch: removes monomorphic variants", { # Store genotype matrix directly in genotype_handle mock_ld_data <- new( "LdData", + variants_gr, correlation = NULL, genotypeHandle = X, - variants = variants_gr, snpIdx = NULL, blockMetadata = blockMetadata ) @@ -3573,6 +3573,99 @@ test_that(".ldFromSketch still errors on a genuinely-absent variant after reconc ) }) +test_that(".ldFromSketch rejects an ldSketch that is not a genotype panel", { + expect_error( + pecotmr:::.ldFromSketch( + "not_a_handle", + c("chr1:100:A:G", "chr1:200:A:G") + ), + "ldSketch must be a genotype panel" + ) +}) + +# A synthetic panel + dosage extractor, so the shape of what .ldFromSketch +# returns is checked without a plink2 fixture (and so without pgenlibr). +# @noRd +.lds_makeHandle <- function(snpN = 6L, nSamples = 30L) { + new( + "GenotypeHandle", + path = "/tmp/sketch.gds", + format = "gds", + snpInfo = data.frame( + SNP = sprintf("chr1:%d:A:G", 100L * seq_len(snpN)), + CHR = rep("1", snpN), + BP = seq(100L, by = 100L, length.out = snpN), + A1 = rep("A", snpN), + A2 = rep("G", snpN), + stringsAsFactors = FALSE + ), + nSamples = nSamples, + sampleIds = sprintf("s%d", seq_len(nSamples)), + pgenPtr = NULL + ) +} + +# @noRd +.lds_mockExtractor <- function(seed = 7, nSamples = 30L) { + function(handle, snpIdx, meanImpute = TRUE) { + set.seed(seed) + nSnp <- nrow(getSnpInfo(handle)) + panel <- matrix( + rbinom(nSamples * nSnp, 2, 0.3), + nrow = nSamples, + ncol = nSnp, + dimnames = list(getSampleIds(handle), getSnpInfo(handle)$SNP) + ) + rr <- GenomicRanges::GRanges( + seqnames = str_c("chr", getSnpInfo(handle)$CHR[snpIdx]), + ranges = IRanges::IRanges( + start = getSnpInfo(handle)$BP[snpIdx], + width = 1L + ) + ) + S4Vectors::mcols(rr) <- S4Vectors::DataFrame( + SNP = getSnpInfo(handle)$SNP[snpIdx], + A1 = getSnpInfo(handle)$A1[snpIdx], + A2 = getSnpInfo(handle)$A2[snpIdx] + ) + dosage <- t(panel[, snpIdx, drop = FALSE]) + dimnames(dosage) <- list( + getSnpInfo(handle)$SNP[snpIdx], + getSampleIds(handle) + ) + SummarizedExperiment::SummarizedExperiment( + assays = list(dosage = dosage), + rowRanges = rr, + colData = S4Vectors::DataFrame( + sampleId = getSampleIds(handle), + row.names = getSampleIds(handle) + ) + ) + } +} + +test_that(".ldFromSketch returns a symmetric unit-diagonal LD matrix", { + h <- .lds_makeHandle() + local_mocked_bindings( + extractBlockGenotypes = .lds_mockExtractor(), + .package = "pecotmr" + ) + ids <- c("chr1:200:A:G", "chr1:400:A:G", "chr1:500:A:G") + R <- pecotmr:::.ldFromSketch(h, ids) + expect_true(is.matrix(R)) + expect_equal(dimnames(R), list(ids, ids)) + expect_equal(unname(diag(R)), rep(1, 3), tolerance = 1e-12) + expect_equal(R, t(R), tolerance = 1e-12) +}) + +test_that(".ldFromSketch errors on a variant the panel does not carry", { + h <- .lds_makeHandle() + expect_error( + pecotmr:::.ldFromSketch(h, c("chr1:100:A:G", "ghost")), + "variant id.*not present in the LD sketch" + ) +}) + # ============================================================================= # Additional coverage: .requireMatchingLdSketches error paths # ============================================================================= @@ -3599,11 +3692,11 @@ test_that(".requireMatchingLdSketches errors when panels differ in a column", { si2$A1[1] <- if (identical(si2$A1[1], "A")) "C" else "A" # mutate one allele h2 <- new( "GenotypeHandle", - path = h@path, - format = h@format, + path = getPath(h), + format = getFormat(h), snpInfo = si2, - nSamples = h@nSamples, - sampleIds = h@sampleIds, + nSamples = getNSamples(h), + sampleIds = getSampleIds(h), pgenPtr = NULL, chromPaths = character(0) ) diff --git a/tests/testthat/test_ldData.R b/tests/testthat/test_ldData.R index 2c3e3593..2642b96c 100644 --- a/tests/testthat/test_ldData.R +++ b/tests/testthat/test_ldData.R @@ -46,28 +46,28 @@ context("LdData accessors") rbinom(n_samples * nrow(getSnpInfo(handle)), 2, 0.3), nrow = n_samples, ncol = nrow(getSnpInfo(handle)), - dimnames = list(handle@sampleIds, getSnpInfo(handle)$SNP) + dimnames = list(getSampleIds(handle), getSnpInfo(handle)$SNP) ) sub <- panel[, snpIdx, drop = FALSE] rr <- GenomicRanges::GRanges( - seqnames = paste0("chr", handle@snpInfo$CHR[snpIdx]), + seqnames = paste0("chr", getSnpInfo(handle)$CHR[snpIdx]), ranges = IRanges::IRanges( - start = handle@snpInfo$BP[snpIdx], + start = getSnpInfo(handle)$BP[snpIdx], width = 1L ) ) S4Vectors::mcols(rr) <- S4Vectors::DataFrame( SNP = getSnpInfo(handle)$SNP[snpIdx], - A1 = handle@snpInfo$A1[snpIdx], - A2 = handle@snpInfo$A2[snpIdx] + A1 = getSnpInfo(handle)$A1[snpIdx], + A2 = getSnpInfo(handle)$A2[snpIdx] ) cd <- S4Vectors::DataFrame( - sampleId = handle@sampleIds, - row.names = handle@sampleIds + sampleId = getSampleIds(handle), + row.names = getSampleIds(handle) ) dosage <- t(sub) rownames(dosage) <- getSnpInfo(handle)$SNP[snpIdx] - colnames(dosage) <- handle@sampleIds + colnames(dosage) <- getSampleIds(handle) SummarizedExperiment::SummarizedExperiment( assays = list(dosage = dosage), rowRanges = rr, @@ -167,10 +167,10 @@ test_that("getCorrelation: mixture handles without mixtureWeights errors", { # Build via new() to skip the constructor's mixtureWeights validity check. ld <- new( "LdData", + .ld_makeVariants(), correlation = NULL, genotypeHandle = list(gh, gh), snpIdx = 1:4, - variants = .ld_makeVariants(), blockMetadata = S4Vectors::DataFrame(x = 1), nRef = 0L, mixtureWeights = NULL @@ -186,10 +186,10 @@ test_that("getCorrelation: mixture panels of differing dim error", { gh <- .ld_makeHandle(snp_n = 4L) ld <- new( "LdData", + .ld_makeVariants(), correlation = NULL, genotypeHandle = list(gh_small, gh), snpIdx = 1:3, - variants = .ld_makeVariants(), blockMetadata = S4Vectors::DataFrame(x = 1), nRef = 0L, mixtureWeights = c(0.5, 0.5) @@ -239,10 +239,10 @@ test_that("getGenotypes: matrix handle is returned unchanged", { ) ld <- new( "LdData", + .ld_makeVariants(), correlation = NULL, genotypeHandle = X, snpIdx = NULL, - variants = .ld_makeVariants(), blockMetadata = S4Vectors::DataFrame(x = 1), nRef = 0L, mixtureWeights = NULL @@ -389,7 +389,7 @@ test_that("LdData validation rejects empty variants", { gr <- GenomicRanges::GRanges() expect_error( LdData(correlation = R, variants = gr, blockMetadata = data.frame()), - "must not be empty" + "must cover >= 1 variant" ) }) @@ -555,7 +555,7 @@ test_that("LdData: validity rejects empty variants", { variants = GenomicRanges::GRanges(), blockMetadata = S4Vectors::DataFrame(x = 1) ), - "'variants' must not be empty" + "must cover >= 1 variant" ) }) @@ -839,8 +839,11 @@ test_that("the payload slots refuse values of the wrong shape", { test_that("blockMetadata takes ranges as well as tables", { gr <- GenomicRanges::GRanges("chr1", IRanges::IRanges(1L, 1000L)) - for (v in list(gr, data.frame(chrom = "chr1"), - S4Vectors::DataFrame(chrom = "chr1"))) { + for (v in list( + gr, + data.frame(chrom = "chr1"), + S4Vectors::DataFrame(chrom = "chr1") + )) { ld <- LdData( correlation = matrix(1), variants = .ld_makeVariants(), diff --git a/tests/testthat/test_manifestLoaders.R b/tests/testthat/test_manifestLoaders.R index 20b04d1b..64fe550b 100644 --- a/tests/testthat/test_manifestLoaders.R +++ b/tests/testthat/test_manifestLoaders.R @@ -122,7 +122,7 @@ test_that(".reconcileScalar resolves arg/column and flags conflicts", { test_that(".detectGenotypeFormat builds a PLINK1 handle from a prefix", { h <- pecotmr:::.detectGenotypeFormat(.toyRefPrefix()) expect_s4_class(h, "GenotypeHandle") - expect_equal(h@format, "plink1") + expect_equal(getFormat(h), "plink1") }) test_that("BCF sumstats are rejected", { @@ -197,7 +197,7 @@ test_that("loadGwasSumStatsFromManifest builds from a data.frame manifest", { c("SNP", "A1", "A2", "Z", "N") %in% colnames(S4Vectors::mcols(obj[[1L]])) )) - expect_equal(length(obj@qcInfo), 0L) # loaders run no QC + expect_equal(length(getQcInfo(obj)), 0L) # loaders run no QC }) test_that("loadGwasSumStatsFromManifest reads a manifest file and reconciles genome", { @@ -662,7 +662,7 @@ test_that("loadMultiStudyQtlDatasetFromManifest builds from >=2 studies", { msd <- loadMultiStudyQtlDatasetFromManifest(manifest) expect_s4_class(msd, "MultiStudyQtlDataset") expect_true(methods::validObject(msd)) - expect_equal(sort(names(msd@qtlDatasets)), c("study1", "study2")) + expect_equal(sort(names(getQtlDatasets(msd))), c("study1", "study2")) }) test_that("loadMultiStudyQtlDatasetFromManifest attaches a summary-only study", { @@ -691,7 +691,7 @@ test_that("loadMultiStudyQtlDatasetFromManifest attaches a summary-only study", ) expect_true(methods::validObject(msd)) expect_s4_class(msd@sumStats, "QtlSumStats") - expect_equal(names(msd@qtlDatasets), "study1") + expect_equal(names(getQtlDatasets(msd)), "study1") }) # =========================================================================== @@ -752,13 +752,13 @@ test_that(".detectGenotypeFormat dispatches by extension", { "Could not determine genotype format" ) h <- pecotmr:::.detectGenotypeFormat(paste0(.toyRefPrefix(), ".bed")) - expect_equal(h@format, "plink1") + expect_equal(getFormat(h), "plink1") }) test_that(".resolveLdSketch accepts a genoMeta vector, a path, and rejects bad input", { sharded <- pecotmr:::.resolveLdSketch(c("22" = .toyRefPrefix())) expect_s4_class(sharded, "GenotypeHandle") - expect_true("22" %in% names(sharded@chromPaths)) + expect_true("22" %in% names(getChromPaths(sharded))) expect_s4_class( pecotmr:::.resolveLdSketch(.toyRefPrefix()), "GenotypeHandle" @@ -807,7 +807,7 @@ test_that(".materializeLdSketch passes a handle through and restricts spec shard spec <- c("22" = .toyRefPrefix(), "21" = "/no/such/chr21/prefix") restricted <- pecotmr:::.materializeLdSketch(spec, "22") expect_s4_class(restricted, "GenotypeHandle") - expect_equal(names(restricted@chromPaths), "22") + expect_equal(names(getChromPaths(restricted)), "22") expect_error(pecotmr:::.materializeLdSketch(spec, c("22", "21"))) }) diff --git a/tests/testthat/test_mashWrapper.R b/tests/testthat/test_mashWrapper.R index 0c52d19f..649cc612 100644 --- a/tests/testthat/test_mashWrapper.R +++ b/tests/testthat/test_mashWrapper.R @@ -1619,8 +1619,8 @@ test_that("qtlSumStatsFromZMatrix: placeholders and pass-through qcInfo are set" expect_equal(unique(mc$A1), "T") expect_equal(unique(mc$A2), "C") expect_equal(unique(mc$N), 500L) - expect_equal(qss@qcInfo$role, "strong") - expect_equal(length(qss@qcInfo$entryAudit), 2L) # one slot per context + expect_equal(getQcInfo(qss)$role, "strong") + expect_equal(length(getQcInfo(qss)$entryAudit), 2L) # one slot per context }) test_that("qtlSumStatsFromZMatrix: columns can map to traits or context x trait pairs", { @@ -1767,7 +1767,7 @@ test_that("qtlSumStatsFromBetaMatrix: NULL rownames -> synthetic ids; placeholde expect_equal(mc$SNP, c("var1", "var2")) expect_equal(unique(mc$A1), "T") expect_equal(unique(mc$N), 500L) - expect_equal(qss@qcInfo$role, "strong") + expect_equal(getQcInfo(qss)$role, "strong") }) # =========================================================================== diff --git a/tests/testthat/test_qtlAssociationPostprocess.R b/tests/testthat/test_qtlAssociationPostprocess.R index ac69d48a..e992c3c4 100644 --- a/tests/testthat/test_qtlAssociationPostprocess.R +++ b/tests/testthat/test_qtlAssociationPostprocess.R @@ -81,7 +81,7 @@ test_that("qtlAssociationPostprocess enriches with package-computed columns", { ) expect_s4_class(r, "QtlSumStats") - expect_false(is.null(r@qcInfo$associationPostprocess)) # recipe stashed + expect_false(is.null(getQcInfo(r)$associationPostprocess)) # recipe stashed # Bonferroni original == min over variants of p.adjust(P, "bonferroni", n). expP <- vapply( diff --git a/tests/testthat/test_qtlEnrichmentPipeline.R b/tests/testthat/test_qtlEnrichmentPipeline.R index f9f3885f..87cbf317 100644 --- a/tests/testthat/test_qtlEnrichmentPipeline.R +++ b/tests/testthat/test_qtlEnrichmentPipeline.R @@ -153,7 +153,7 @@ test_that("qtlEnrichmentPipeline: ldSketch mismatch errors", { # Build the QTL with a sketch carrying a different sample set. gfmr <- .qep_makeGwasFmr() qSketch <- .qep_makeHandle() - qSketch@sampleIds <- paste0("z", seq_len(qSketch@nSamples)) + qSketch@sampleIds <- paste0("z", seq_len(getNSamples(qSketch))) qfmr <- QtlFineMappingResult( study = "Q1", context = "c1", diff --git a/tests/testthat/test_sldscWrapper.R b/tests/testthat/test_sldscWrapper.R index 699380b6..bdf4bbb9 100644 --- a/tests/testthat/test_sldscWrapper.R +++ b/tests/testthat/test_sldscWrapper.R @@ -9,7 +9,6 @@ # (computeSldscMRef/computeSldscAnnotSd/isBinarySldscAnnot/standardizeSldscTrait) # operate on an in-memory SldscData and are tested with in-memory fixtures. - # ============================================================================= # .sldscDetectAnnotCols # ============================================================================= diff --git a/tests/testthat/test_sumstatsQc.R b/tests/testthat/test_sumstatsQc.R index 2098b455..88166b31 100644 --- a/tests/testthat/test_sumstatsQc.R +++ b/tests/testthat/test_sumstatsQc.R @@ -3084,31 +3084,31 @@ context("summaryStatsQc") function(handle, snpIdx, meanImpute = TRUE) { set.seed(seed) panel <- matrix( - rbinom(n_samples * nrow(handle@snpInfo), 2, 0.3), + rbinom(n_samples * nrow(getSnpInfo(handle)), 2, 0.3), nrow = n_samples, ncol = nrow(getSnpInfo(handle)), - dimnames = list(handle@sampleIds, getSnpInfo(handle)$SNP) + dimnames = list(getSampleIds(handle), getSnpInfo(handle)$SNP) ) sub <- panel[, snpIdx, drop = FALSE] rr <- GenomicRanges::GRanges( - seqnames = paste0("chr", handle@snpInfo$CHR[snpIdx]), + seqnames = paste0("chr", getSnpInfo(handle)$CHR[snpIdx]), ranges = IRanges::IRanges( - start = handle@snpInfo$BP[snpIdx], + start = getSnpInfo(handle)$BP[snpIdx], width = 1L ) ) S4Vectors::mcols(rr) <- S4Vectors::DataFrame( - SNP = handle@snpInfo$SNP[snpIdx], - A1 = handle@snpInfo$A1[snpIdx], - A2 = handle@snpInfo$A2[snpIdx] + SNP = getSnpInfo(handle)$SNP[snpIdx], + A1 = getSnpInfo(handle)$A1[snpIdx], + A2 = getSnpInfo(handle)$A2[snpIdx] ) cd <- S4Vectors::DataFrame( - sampleId = handle@sampleIds, - row.names = handle@sampleIds + sampleId = getSampleIds(handle), + row.names = getSampleIds(handle) ) dosage <- t(sub) - rownames(dosage) <- handle@snpInfo$SNP[snpIdx] - colnames(dosage) <- handle@sampleIds + rownames(dosage) <- getSnpInfo(handle)$SNP[snpIdx] + colnames(dosage) <- getSampleIds(handle) SummarizedExperiment::SummarizedExperiment( assays = list(dosage = dosage), rowRanges = rr, @@ -3209,7 +3209,7 @@ test_that("summaryStatsQc: PIP screen off leaves the harmonized set intact", { # as real genotype LD sketches do, so a re-keyed entry SNP resolves against it. .ssQ_makeHandleVid <- function(snp_n = 8L, n_samples = 60L) { h <- .ssQ_makeHandle(snp_n, n_samples) - h@snpInfo$SNP <- paste0("chr1:", h@snpInfo$BP, ":G:A") + h@snpInfo$SNP <- paste0("chr1:", getSnpInfo(h)$BP, ":G:A") h } @@ -3297,7 +3297,7 @@ test_that("summaryStatsQc: zMismatchQc reconciles a chr-prefix difference vs the # canonical chr-prefixed form, so the opt-in z-mismatch panel match must # reconcile the prefix (previously errored "absent from the ldSketch panel"). h <- .ssQ_makeHandle() - h@snpInfo$SNP <- paste0("1:", h@snpInfo$BP, ":G:A") # non-chr-prefixed + h@snpInfo$SNP <- paste0("1:", getSnpInfo(h)$BP, ":G:A") # non-chr-prefixed ss <- GwasSumStats( study = "g1", entry = list(.ssQ_makeEntryGr()), @@ -3790,7 +3790,13 @@ test_that("summaryStatsQc: impute scopes the reference panel/dosage to the regio # dosage must be materialized for just those 4 panel variants -- NOT the whole # 8-variant sketch (the bug that makes --impute unusable on a per-chromosome # sketch: it built dosage for seq_len(nrow(sketchSnpInfo))). + # + # Observed at extractBlockGenotypes: the single reader every dosage path + # funnels through, including the DelayedArray seed behind assay(). The + # panel reads its own dosage through the assay now, so watching + # .dosageMatrix would no longer see the RAISS read at all. cap <- new.env(parent = emptyenv()) + cap$idx <- list() ss <- GwasSumStats( study = "g1", entry = list(.ssQ_makeEntryGr( @@ -3800,17 +3806,20 @@ test_that("summaryStatsQc: impute scopes the reference panel/dosage to the regio genome = "hg19", ldSketch = .ssQ_makeHandle(snp_n = 8L, n_samples = 60L) ) + inner <- .ssQ_mockExtractor(n_samples = 60L) local_mocked_bindings( - .dosageMatrix = function(handle, snpIdx, meanImpute = TRUE) { - cap$snpIdx <- snpIdx - matrix(0, nrow = handle@nSamples, ncol = length(snpIdx)) + extractBlockGenotypes = function(handle, snpIdx, meanImpute = TRUE) { + cap$idx <- c(cap$idx, list(snpIdx)) + inner(handle, snpIdx, meanImpute) }, raiss = function(...) NULL, .package = "pecotmr" ) suppressWarnings(summaryStatsQc(ss, impute = TRUE, nCutoff = 0)) - expect_equal(length(cap$snpIdx), 4L) # region window, not the full 8-SNP sketch - expect_true(all(cap$snpIdx %in% 1:4)) # only the in-window panel indices + expect_gt(length(cap$idx), 0L) + # No read reaches past the region window's 4 panel variants. + expect_true(all(map_lgl(cap$idx, function(i) all(i %in% 1:4)))) + expect_equal(max(map_int(cap$idx, length)), 4L) }) test_that("summaryStatsQc: impute = TRUE with raiss returning NULL records 0 imputed", { @@ -3994,31 +4003,31 @@ context("sumstatsQc internal helpers") function(handle, snpIdx, meanImpute = TRUE) { set.seed(seed) panel <- matrix( - rbinom(n_samples * nrow(handle@snpInfo), 2, 0.3), + rbinom(n_samples * nrow(getSnpInfo(handle)), 2, 0.3), nrow = n_samples, ncol = nrow(getSnpInfo(handle)), - dimnames = list(handle@sampleIds, getSnpInfo(handle)$SNP) + dimnames = list(getSampleIds(handle), getSnpInfo(handle)$SNP) ) sub <- panel[, snpIdx, drop = FALSE] rr <- GenomicRanges::GRanges( - seqnames = paste0("chr", handle@snpInfo$CHR[snpIdx]), + seqnames = paste0("chr", getSnpInfo(handle)$CHR[snpIdx]), ranges = IRanges::IRanges( - start = handle@snpInfo$BP[snpIdx], + start = getSnpInfo(handle)$BP[snpIdx], width = 1L ) ) S4Vectors::mcols(rr) <- S4Vectors::DataFrame( - SNP = handle@snpInfo$SNP[snpIdx], - A1 = handle@snpInfo$A1[snpIdx], - A2 = handle@snpInfo$A2[snpIdx] + SNP = getSnpInfo(handle)$SNP[snpIdx], + A1 = getSnpInfo(handle)$A1[snpIdx], + A2 = getSnpInfo(handle)$A2[snpIdx] ) cd <- S4Vectors::DataFrame( - sampleId = handle@sampleIds, - row.names = handle@sampleIds + sampleId = getSampleIds(handle), + row.names = getSampleIds(handle) ) dosage <- t(sub) - rownames(dosage) <- handle@snpInfo$SNP[snpIdx] - colnames(dosage) <- handle@sampleIds + rownames(dosage) <- getSnpInfo(handle)$SNP[snpIdx] + colnames(dosage) <- getSampleIds(handle) SummarizedExperiment::SummarizedExperiment( assays = list(dosage = dosage), rowRanges = rr, @@ -4406,6 +4415,119 @@ test_that(".applyLdMismatchQcToEntry: errors on variants absent from the sketch" ) }) +# =========================================================================== +# The RSS QC steps agree with the panel about which variants exist +# +# Harmonization decides the surviving variant set; kriging, LD-mismatch QC and +# the fine-mapping LD build then assert that every survivor is in the panel. +# A panel carrying a variant AND its own flip used to break that: harmonization +# emitted the variant twice (opposite Z, one per orientation) and the lookups, +# which go through matchVariants, refused both and aborted. +# =========================================================================== + +# @noRd +.ssTwin_makeHandle <- function(nSamples = 60L) { + new( + "GenotypeHandle", + path = "/tmp/sketch.gds", + format = "gds", + snpInfo = data.frame( + SNP = c( + "chr1:100:A:G", + "chr1:200:A:G", + "chr1:200:G:A", + "chr1:300:C:T" + ), + CHR = rep("1", 4L), + BP = c(100L, 200L, 200L, 300L), + A1 = c("G", "G", "A", "T"), + A2 = c("A", "A", "G", "C"), + stringsAsFactors = FALSE + ), + nSamples = nSamples, + sampleIds = paste0("s", seq_len(nSamples)), + pgenPtr = NULL + ) +} + +# @noRd +.ssTwin_df <- function() { + data.frame( + chrom = rep("1", 3L), + pos = c(100L, 200L, 300L), + SNP = c("chr1:100:A:G", "chr1:200:A:G", "chr1:300:C:T"), + A1 = c("G", "G", "T"), + A2 = c("A", "A", "C"), + Z = c(1.1, 2.2, 3.3), + N = rep(1000, 3L), + stringsAsFactors = FALSE + ) +} + +# @noRd +.ssTwin_opts <- function() { + list( + matchMinProp = 0, + removeIndels = FALSE, + removeStrandAmbiguous = TRUE, + alleleFlipKriging = TRUE, + nForPip = 1000, + zMismatchQc = "slalom" + ) +} + +test_that("harmonization drops a twin-pair variant instead of duplicating it", { + harm <- suppressMessages(pecotmr:::.qcHarmonizeEntry( + .ssTwin_df(), + .ssTwin_makeHandle(), + .ssTwin_opts(), + NA_character_ + )) + # 3 in, 2 out -- never 4: the undecidable chr1:200 leaves, it does not + # come back twice with opposite signs. + expect_equal(nrow(harm$df), 2L) + expect_false(any(harm$df$pos == 200L)) + expect_equal(harm$counts$harmDropped, 1L) +}) + +test_that("kriging and LD-mismatch QC accept every harmonized variant", { + handle <- .ssTwin_makeHandle() + opts <- .ssTwin_opts() + local_mocked_bindings( + extractBlockGenotypes = .ssQ_mockExtractor(n_samples = 60L), + .package = "pecotmr" + ) + harm <- suppressMessages( + pecotmr:::.qcHarmonizeEntry(.ssTwin_df(), handle, opts, NA_character_) + ) + expect_no_error(suppressMessages( + pecotmr:::.qcKrigingFlip(harm$df, handle, opts, NA_character_) + )) + expect_no_error(suppressMessages( + pecotmr:::.qcMismatchQc(harm$df, handle, opts, NA_character_) + )) +}) + +test_that("the fine-mapping LD build accepts every harmonized variant", { + handle <- .ssTwin_makeHandle() + local_mocked_bindings( + extractBlockGenotypes = .ssQ_mockExtractor(n_samples = 60L), + .package = "pecotmr" + ) + harm <- suppressMessages(pecotmr:::.qcHarmonizeEntry( + .ssTwin_df(), + handle, + .ssTwin_opts(), + NA_character_ + )) + R <- pecotmr:::.ldFromSketch( + handle, + as.character(harm$df$SNP), + label = "fineMappingPipeline" + ) + expect_equal(dim(R), c(2L, 2L)) +}) + test_that(".applyLdMismatchQcToEntry: NA outlier flags from slalom are kept (not dropped)", { # Regression test: slalom (and dentist on degenerate inputs) can leave # NA in the `outlier` column for variants whose per-variant statistic @@ -7191,33 +7313,33 @@ test_that("summaryStatsQc kriging QC sign-flips an LD-inconsistent variant and r # the row. corrExtractor <- function(handle, snpIdx, meanImpute = TRUE) { set.seed(42) - n <- length(handle@sampleIds) + n <- length(getSampleIds(handle)) k <- length(snpIdx) f <- rnorm(n) # shared latent factor M <- sapply(seq_len(k), function(j) { sqrt(0.7) * f + sqrt(0.3) * rnorm(n) }) rr <- GenomicRanges::GRanges( - seqnames = paste0("chr", handle@snpInfo$CHR[snpIdx]), + seqnames = paste0("chr", getSnpInfo(handle)$CHR[snpIdx]), ranges = IRanges::IRanges( - start = handle@snpInfo$BP[snpIdx], + start = getSnpInfo(handle)$BP[snpIdx], width = 1L ) ) S4Vectors::mcols(rr) <- S4Vectors::DataFrame( - SNP = handle@snpInfo$SNP[snpIdx], - A1 = handle@snpInfo$A1[snpIdx], - A2 = handle@snpInfo$A2[snpIdx] + SNP = getSnpInfo(handle)$SNP[snpIdx], + A1 = getSnpInfo(handle)$A1[snpIdx], + A2 = getSnpInfo(handle)$A2[snpIdx] ) dosage <- t(M) - rownames(dosage) <- handle@snpInfo$SNP[snpIdx] - colnames(dosage) <- handle@sampleIds + rownames(dosage) <- getSnpInfo(handle)$SNP[snpIdx] + colnames(dosage) <- getSampleIds(handle) SummarizedExperiment::SummarizedExperiment( assays = list(dosage = dosage), rowRanges = rr, colData = S4Vectors::DataFrame( - sampleId = handle@sampleIds, - row.names = handle@sampleIds + sampleId = getSampleIds(handle), + row.names = getSampleIds(handle) ) ) } @@ -7729,6 +7851,46 @@ test_that(".raissFlipPairMask needs both orientations, not just an indel", { expect_equal(.raissFlipPairMask(ids), c(TRUE, TRUE, FALSE, FALSE)) }) +test_that(".raissFlipPairMask sees a strand-recorded second orientation", { + # A:G and C:T are one variant written on opposite strands AND swapped, so + # matchVariants() refuses the position. A literal A1/A2 reversal key read + # this as two unrelated variants and let RAISS impute the pair back. + ids <- normalizeVariantId(c("chr1:100:A:G", "chr1:100:C:T")) + expect_equal(.raissFlipPairMask(ids), c(TRUE, TRUE)) +}) + +test_that(".raissFlipPairMask agrees with matchVariants on what is ambiguous", { + # The mask exists to stop imputation putting back what matchVariants drops, + # so the two have to answer the same question the same way. + ids <- normalizeVariantId(c( + "chr1:100:A:G", + "chr1:100:C:T", + "chr1:200:A:G", + "chr1:200:G:A", + "chr1:300:A:G", + "chr1:400:A:T" + )) + m <- matchVariants(ids, ids, removeStrandAmbiguous = FALSE) + expect_equal(.raissFlipPairMask(ids), !is_in(seq_along(ids), m$idxA)) + expect_equal( + .raissFlipPairMask(ids), + c(TRUE, TRUE, TRUE, TRUE, FALSE, FALSE) + ) +}) + +test_that(".raissFlipPairMask spares the first of two identical entries", { + # Nothing is ambiguous about a variant listed twice -- the first copy is + # imputable. The second is not: imputing it would put one id into the + # sumstats twice, and matchVariants answers a repeated id once, so the + # next LD lookup would abort on the copy it could not place. + ids <- normalizeVariantId(c("chr1:100:A:G", "chr1:100:A:G")) + expect_equal(.raissFlipPairMask(ids), c(FALSE, TRUE)) +}) + +test_that(".raissFlipPairMask handles an empty panel", { + expect_equal(.raissFlipPairMask(character(0)), logical(0)) +}) + # =========================================================================== # RAISS: positions the GWAS already typed (ported from PR #573) # diff --git a/tests/testthat/test_twasWeightsPipeline.R b/tests/testthat/test_twasWeightsPipeline.R index bebb87f4..17b8c9a6 100644 --- a/tests/testthat/test_twasWeightsPipeline.R +++ b/tests/testthat/test_twasWeightsPipeline.R @@ -95,31 +95,31 @@ context("twasWeightsPipeline (S4 dispatch) with mocked weight methods") function(handle, snpIdx, meanImpute = TRUE) { set.seed(seed) panel <- matrix( - rbinom(n_samples * nrow(handle@snpInfo), 2, 0.3), + rbinom(n_samples * nrow(getSnpInfo(handle)), 2, 0.3), nrow = n_samples, - ncol = nrow(handle@snpInfo), - dimnames = list(handle@sampleIds, handle@snpInfo$SNP) + ncol = nrow(getSnpInfo(handle)), + dimnames = list(getSampleIds(handle), getSnpInfo(handle)$SNP) ) sub <- panel[, snpIdx, drop = FALSE] rr <- GenomicRanges::GRanges( - seqnames = paste0("chr", handle@snpInfo$CHR[snpIdx]), + seqnames = paste0("chr", getSnpInfo(handle)$CHR[snpIdx]), ranges = IRanges::IRanges( - start = handle@snpInfo$BP[snpIdx], + start = getSnpInfo(handle)$BP[snpIdx], width = 1L ) ) S4Vectors::mcols(rr) <- S4Vectors::DataFrame( - SNP = handle@snpInfo$SNP[snpIdx], - A1 = handle@snpInfo$A1[snpIdx], - A2 = handle@snpInfo$A2[snpIdx] + SNP = getSnpInfo(handle)$SNP[snpIdx], + A1 = getSnpInfo(handle)$A1[snpIdx], + A2 = getSnpInfo(handle)$A2[snpIdx] ) cd <- S4Vectors::DataFrame( - sampleId = handle@sampleIds, - row.names = handle@sampleIds + sampleId = getSampleIds(handle), + row.names = getSampleIds(handle) ) dosage <- t(sub) - rownames(dosage) <- handle@snpInfo$SNP[snpIdx] - colnames(dosage) <- handle@sampleIds + rownames(dosage) <- getSnpInfo(handle)$SNP[snpIdx] + colnames(dosage) <- getSampleIds(handle) SummarizedExperiment::SummarizedExperiment( assays = list(dosage = dosage), rowRanges = rr, @@ -2149,101 +2149,6 @@ test_that(".twasFineMappingFits: ignores non-susie methods (e.g. lasso)", { expect_equal(names(out), "susie") }) -# =========================================================================== -# .twasLdFromSketch (mocked extractBlockGenotypes) -# =========================================================================== - -.tw_makeSketchHandle <- function(snp_n = 6L, n_samples = 30L) { - new( - "GenotypeHandle", - path = "/tmp/sketch.gds", - format = "gds", - snpInfo = data.frame( - SNP = sprintf("chr1:%d:A:G", 100L * (seq_len(snp_n))), - CHR = rep("1", snp_n), - BP = seq(100L, by = 100L, length.out = snp_n), - A1 = rep("A", snp_n), - A2 = rep("G", snp_n), - stringsAsFactors = FALSE - ), - nSamples = n_samples, - sampleIds = paste0("s", seq_len(n_samples)), - pgenPtr = NULL - ) -} - -.tw_mockExtractor <- function(seed = 7, n_samples = 30L) { - function(handle, snpIdx, meanImpute = TRUE) { - set.seed(seed) - panel <- matrix( - rbinom(n_samples * nrow(handle@snpInfo), 2, 0.3), - nrow = n_samples, - ncol = nrow(handle@snpInfo), - dimnames = list(handle@sampleIds, handle@snpInfo$SNP) - ) - sub <- panel[, snpIdx, drop = FALSE] - rr <- GenomicRanges::GRanges( - seqnames = paste0("chr", handle@snpInfo$CHR[snpIdx]), - ranges = IRanges::IRanges( - start = handle@snpInfo$BP[snpIdx], - width = 1L - ) - ) - S4Vectors::mcols(rr) <- S4Vectors::DataFrame( - SNP = handle@snpInfo$SNP[snpIdx], - A1 = handle@snpInfo$A1[snpIdx], - A2 = handle@snpInfo$A2[snpIdx] - ) - cd <- S4Vectors::DataFrame( - sampleId = handle@sampleIds, - row.names = handle@sampleIds - ) - dosage <- t(sub) - rownames(dosage) <- handle@snpInfo$SNP[snpIdx] - colnames(dosage) <- handle@sampleIds - SummarizedExperiment::SummarizedExperiment( - assays = list(dosage = dosage), - rowRanges = rr, - colData = cd - ) - } -} - -test_that(".twasLdFromSketch: rejects a non-panel ldSketch", { - expect_error( - pecotmr:::.twasLdFromSketch( - "not_a_handle", - c("chr1:100:A:G", "chr1:200:A:G") - ), - "ldSketch must be a genotype panel" - ) -}) - -test_that(".twasLdFromSketch: variants not in the panel error", { - h <- .tw_makeSketchHandle() - expect_error( - pecotmr:::.twasLdFromSketch(h, c("chr1:100:A:G", "ghost")), - "variant id.*not present in the LD sketch" - ) -}) - -test_that(".twasLdFromSketch: returns a square LD matrix named by variantIds", { - h <- .tw_makeSketchHandle() - local_mocked_bindings( - extractBlockGenotypes = .tw_mockExtractor(), - .package = "pecotmr" - ) - ids <- c("chr1:200:A:G", "chr1:400:A:G", "chr1:500:A:G") - R <- pecotmr:::.twasLdFromSketch(h, ids) - expect_true(is.matrix(R)) - expect_equal(dim(R), c(3L, 3L)) - expect_equal(rownames(R), ids) - expect_equal(colnames(R), ids) - # Symmetric, diagonal == 1 (sample correlation). - expect_equal(unname(diag(R)), c(1, 1, 1), tolerance = 1e-12) - expect_equal(R, t(R), tolerance = 1e-12) -}) - # =========================================================================== # .twasWeightsPipelineMatrix: susieFit pre-fit pass-through # =========================================================================== @@ -3516,7 +3421,7 @@ test_that(".solveEnsembleGlmnet: solver failure falls back to equal weights", { # @noRd .twrssf_qcd <- function() { data(qtlSumStatsExample, envir = environment()) - suppressWarnings(suppressMessages(summaryStatsQc(qtlSumStatsExample))) + suppressMessages(summaryStatsQc(qtlSumStatsExample)) } # @noRd diff --git a/tests/testthat/test_variantId.R b/tests/testthat/test_variantId.R index 3cb30c2c..14994588 100644 --- a/tests/testthat/test_variantId.R +++ b/tests/testthat/test_variantId.R @@ -877,10 +877,14 @@ test_that("matchVariants allowFlip = FALSE matches exact alleles only (no swap)" } test_that("harmonizeAlleles warns and returns empty when nothing overlaps", { - res <- suppressWarnings(pecotmr:::harmonizeAlleles( - .vid_df("1", 100, "A", "G"), - .vid_df("1", 999, "A", "G") - )) + # Asserted, not suppressed: the warning is half of what this test names. + expect_warning( + res <- pecotmr:::harmonizeAlleles( + .vid_df("1", 100, "A", "G"), + .vid_df("1", 999, "A", "G") + ), + "No matching variants found" + ) expect_equal(nrow(res$harmonizedData), 0L) expect_equal(attr(res, "qcCounts")$considered, 0L) }) @@ -991,6 +995,88 @@ test_that("harmonizeAlleles removeIndels = TRUE drops indels", { expect_equal(nrow(res$harmonizedData), 1L) }) +# =========================================================================== +# harmonizeAlleles: one row per target variant +# +# The join is on (chrom, pos), so a reference carrying both orientations of a +# variant used to return that target TWICE with opposite signs -- a phantom +# sign-flipped copy that every matchVariants-based LD lookup downstream then +# refused to match. The resolution now happens here, where the surviving +# variant set is decided. +# =========================================================================== + +test_that("harmonizeAlleles drops a target the reference cannot decide", { + tgt <- .vid_df("1", 100, "A", "G", Z = 2.2) + ref <- .vid_df(c("1", "1"), c(100, 100), c("A", "G"), c("G", "A")) + res <- pecotmr:::harmonizeAlleles( + tgt, + ref, + colToFlip = "Z", + matchMinProp = 0 + ) + expect_equal(nrow(res$harmonizedData), 0L) +}) + +test_that("harmonizeAlleles sees a strand-recorded second orientation too", { + # A:G and C:T are the same ambiguity as A:G and G:A -- only a literal + # string key would miss it. + tgt <- .vid_df("1", 100, "A", "G", Z = 2.2) + ref <- .vid_df(c("1", "1"), c(100, 100), c("A", "C"), c("G", "T")) + res <- pecotmr:::harmonizeAlleles( + tgt, + ref, + colToFlip = "Z", + matchMinProp = 0 + ) + expect_equal(nrow(res$harmonizedData), 0L) +}) + +test_that("harmonizeAlleles keeps one row when the matches agree", { + # A second reference row at the same position with reconcilable alleles + # and the SAME sign is not ambiguous, so the target survives -- once. + tgt <- .vid_df("1", 100, "A", "G", Z = 2.2) + ref <- .vid_df(c("1", "1"), c(100, 100), c("A", "T"), c("G", "C")) + res <- pecotmr:::harmonizeAlleles( + tgt, + ref, + colToFlip = "Z", + matchMinProp = 0 + ) + expect_equal(nrow(res$harmonizedData), 1L) + expect_equal(res$harmonizedData$Z, 2.2) +}) + +test_that("harmonizeAlleles leaves a plain multi-allelic position alone", { + # A:G and A:T at one position are two different variants, not two records + # of one, so the target matches exactly one of them. + tgt <- .vid_df("1", 100, "A", "G", Z = 2.2) + ref <- .vid_df(c("1", "1"), c(100, 100), c("A", "A"), c("G", "T")) + res <- pecotmr:::harmonizeAlleles( + tgt, + ref, + colToFlip = "Z", + matchMinProp = 0 + ) + expect_equal(nrow(res$harmonizedData), 1L) + expect_equal(res$harmonizedData$Z, 2.2) +}) + +test_that("harmonizeAlleles counts what it kept, not what it joined", { + # The counts feed the summaryStatsQc log line; they are taken after the + # resolution, so a dropped target is not also reported as a correction. + tgt <- .vid_df("1", 100, "A", "G", Z = 2.2) + ref <- .vid_df(c("1", "1"), c(100, 100), c("A", "G"), c("G", "A")) + res <- pecotmr:::harmonizeAlleles( + tgt, + ref, + colToFlip = "Z", + matchMinProp = 0 + ) + counts <- attr(res, "qcCounts") + expect_equal(counts$kept, 0L) + expect_equal(counts$signFlip, 0L) +}) + # =========================================================================== # withChrPrefix + asGranges error paths # =========================================================================== @@ -1060,11 +1146,13 @@ test_that("the drop does not depend on the reference's row order", { test_that("a single orientation still matches, flipped or not", { expect_equal( - matchVariants("chr1:100:A:AT", "chr1:100:A:AT")$sign, 1 + matchVariants("chr1:100:A:AT", "chr1:100:A:AT")$sign, + 1 ) # An unambiguous flip is a match, not an error: only disagreement is. expect_equal( - matchVariants("chr1:100:A:AT", "chr1:100:AT:A")$sign, -1 + matchVariants("chr1:100:A:AT", "chr1:100:AT:A")$sign, + -1 ) }) @@ -1079,7 +1167,8 @@ test_that("duplicate reference entries that agree are still matched", { test_that("only the ambiguous target is dropped from a mixed set", { sumstats <- c("chr1:100:A:AT", "chr1:400:A:G", "chr1:500:C:T") panel <- c( - "chr1:100:A:AT", "chr1:100:AT:A", + "chr1:100:A:AT", + "chr1:100:AT:A", "chr1:400:A:G", "chr1:500:T:C" )