diff --git a/.github/environment/pixi.toml b/.github/environment/pixi.toml index 3bb6da8c..45e30589 100644 --- a/.github/environment/pixi.toml +++ b/.github/environment/pixi.toml @@ -9,7 +9,7 @@ platforms = [ [tasks] devtools_document = "cd $GITHUB_WORKSPACE; R -e 'devtools::document()'" -devtools_test = "cd $GITHUB_WORKSPACE; R -e 'devtools::test()'" +devtools_test = "cd $GITHUB_WORKSPACE; R -e 'devtools::test(stop_on_failure = TRUE)'" codecov = "cd $GITHUB_WORKSPACE; R -e 'covr::codecov(quiet = FALSE)'" rcmdcheck = "cd $GITHUB_WORKSPACE; R -e 'rcmdcheck::rcmdcheck()'" bioccheck_git_clone = "cd $GITHUB_WORKSPACE; R -e 'BiocCheck::BiocCheckGitClone()'" diff --git a/R/AllGenerics.R b/R/AllGenerics.R index 039b0f03..1b2f7233 100644 --- a/R/AllGenerics.R +++ b/R/AllGenerics.R @@ -2281,8 +2281,9 @@ setGeneric("credibleSetSummary", function(x, ...) { #' the fitted effect curve and its credible band over the functional grid, one #' row per (credible set, grid position). Wraps the upstream fSuSiE band #' computation (fsusieR's internal \code{update_cal_credible_band.susiF} + -#' \code{get_fitted_effect}). Requires an UNtrimmed fit (the wavelet slots are -#' dropped by trimming); degrades to zero rows for non-fSuSiE or trimmed fits. +#' \code{get_fitted_effect}). Uses saved bands in trimmed fits, or computes +#' them from an untrimmed fit. Returns zero rows when the functional fields +#' are absent. #' #' @param x A \code{FineMappingRow} or \code{FineMappingResultBase}. #' @param ... Ignored. @@ -2303,8 +2304,8 @@ setGeneric("fsusieCredibleBand", function(x, ...) { #' zero (the genomic footprint of each fSuSiE effect), from the upstream #' \code{fsusieR::affected_reg}. Returns a \code{GRanges} with the credible-set #' label, its purity, and the effect \code{direction} -#' (\code{"pos"}/\code{"neg"}; upstream discards it). Requires an UNtrimmed fit; -#' empty otherwise. +#' (\code{"pos"}/\code{"neg"}; upstream discards it). Uses retained functional +#' summaries or an untrimmed fit; empty when those fields are absent. #' #' @param x A \code{FineMappingRow} or \code{FineMappingResultBase}. #' @param ... Ignored. diff --git a/R/GenotypeHandle.R b/R/GenotypeHandle.R index 1b837321..c424eb36 100644 --- a/R/GenotypeHandle.R +++ b/R/GenotypeHandle.R @@ -161,6 +161,7 @@ setMethod("show", "GenotypeHandle", function(object) { #' \code{VariantAnnotation::readVcf} when the source is a VCF. #' @return A \code{GenotypeHandle} object. #' @keywords internal +#' @rdname GenotypeHandle-class GenotypeHandle <- function( path = NULL, plink1Prefix = NULL, diff --git a/R/fineMappingRow.R b/R/fineMappingRow.R index f9ca2bb4..bc2abdd9 100644 --- a/R/fineMappingRow.R +++ b/R/fineMappingRow.R @@ -201,10 +201,10 @@ fineMappingRow <- function(variantIds, susieFit, topLoci, cvResult = NULL) { .isFsusieFit <- function(fit) { !is.null(fit) && is.list(fit) && - !is.null(fit$fitted_wc2) && !is.null(fit$fitted_func) && !is.null(fit$outing_grid) && - !is.null(fit$alpha) + (!is.null(fit$cred_band) || + (!is.null(fit$fitted_wc2) && !is.null(fit$alpha))) } # Chromosome of the fit's variants (matches topLoci's `chrom`, no "chr"). diff --git a/R/fineMappingWrappers.R b/R/fineMappingWrappers.R index 213ba716..5fef0a8e 100644 --- a/R/fineMappingWrappers.R +++ b/R/fineMappingWrappers.R @@ -2238,15 +2238,15 @@ trimFinemappingFit <- function(fit, effectIdx, method, csTables) { } else { common } - # fSuSiE: keep the precomputed variants x features TWAS weight matrix - # (fsusieWeights output, attached as $coef before trimming) so downstream - # TWAS can read it without the dropped wavelet slots. - withCoef <- if (method == "fsusie" && !is.null(fit$coef)) { - list_assign(trimmed, coef = fit$coef) + # fSuSiE: retain the functional summaries alongside the precomputed + # weights, so a trimmed fit still yields curves, bands and TWAS weights + # once the wavelet reconstruction fields are gone. + withFsusie <- if (method == "fsusie") { + .trimAddFsusie(trimmed, fit) } else { trimmed } - `class<-`(withCoef, unique(c(method, "susie"))) + `class<-`(withFsusie, unique(c(method, "susie"))) } # The minimal always-kept subset of a susie fit (pip, credible sets, effect @@ -2341,6 +2341,30 @@ trimFinemappingFit <- function(fit, effectIdx, method, csTables) { # mvsusie-specific slots: per-effect mu2_diag, the coefficient matrix, and the # conditional lfsr array. +# fSuSiE: the fields the wavelet reconstruction would otherwise take with +# it. The base trimmed fit already keeps alpha, lbf_variable and the sets for +# coloc; these carry the fitted curve, its credible band and the functional +# grid. The band is populated BEFORE trimming because +# .fsusiePopulateCredibleBand() needs the reconstruction fields that trimming +# removes. intersect() rather than a per-field guard: an absent field is +# simply not retained. +# @noRd +.trimAddFsusie <- function(trimmed, fit) { + populated <- if (.isFsusieFit(fit)) { + .fsusiePopulateCredibleBand(fit) + } else { + fit + } + keep <- intersect( + c( + "coef", "fitted_func", "cred_band", "outing_grid", "cs", + "csd_X", "lfsr_func" + ), + names(populated) + ) + list_assign(trimmed, !!!populated[keep]) +} + .trimAddMvsusie <- function(trimmed, fit, effectIdx) { list_assign( trimmed, @@ -5550,6 +5574,12 @@ krigingOutlierQc <- function( # @noRd #' @importFrom rlang try_fetch .fsusiePopulateCredibleBand <- function(fit) { + # A native post_processing = "none" run can leave zero-filled placeholder + # bands, so only reuse a band once trimming has removed the + # reconstruction fields -- otherwise recompute it. + if (is.null(fit$fitted_wc2) && !is.null(fit$cred_band)) { + return(fit) + } fn <- try_fetch( get("update_cal_credible_band.susiF", envir = asNamespace("fsusieR")), error = function(cnd) NULL diff --git a/R/jointEngine.R b/R/jointEngine.R index ad64cb90..faa6337e 100644 --- a/R/jointEngine.R +++ b/R/jointEngine.R @@ -277,7 +277,7 @@ setMethod( ) # fsusie joint fit (functional SuSiE over the trait domain; individual-level, -# cross-trait). One per-condition entry per trait, with an optional CV slice. +# cross-trait). Keep the first post-processed entry for this joint fit. # @noRd .jointFitFsusie <- function(group, Xc, Yc, nCond, cfg, args) { if (length(.jgTraitPos(group)) != nCond) { @@ -295,7 +295,7 @@ setMethod( ) raw <- exec(fitFsusie, !!!.splitMethodArgs(fitFsusie, fitArgs)) # Collapse the functional fit to a variants x features weight matrix now - # (trimming later drops fitted_wc/csd_X); store on $coef so a trimmed fit + # (trimming later drops fitted_wc); store on $coef so a trimmed fit # can still yield TWAS weights. fit <- list_assign( raw, @@ -306,14 +306,27 @@ setMethod( ) |> .setFinemappingFitClass("fsusie") cvM <- .jointFsusieCv(Xc, Yc, group, cfg, args, verbose) - map( - seq_len(nCond), - .jointFsusieEntry, + entry <- .jointFsusieEntry( + r = 1L, fit = fit, cvM = cvM, Xc = Xc, cfg = cfg ) + # ONE entry for the whole joint fit, unlike mvsusie's one-per-condition. + # mvsusie's posteriors carry a conditions axis, so .fmSliceCond() gives + # each condition its own 2-D slice and N rows partition the fit at no + # cost. An fsusie fit has no such axis: each trait is a POSITION on a + # shared functional grid, so `fitted_func` / `cred_band` / `outing_grid` + # describe all traits at once and cannot be sliced per trait. One row per + # trait would therefore store N copies of the same curve. + # + # The trait list itself is not lost -- every row carries the full joint + # membership in its `jointTraits` column. Only the grid POSITIONS are + # unrecoverable from that, and they are what maps a curve position back + # to a trait, so they ride on the fit. + entry@susieFit$trait_positions <- .jgTraitPos(group) + list(entry) } # Per-fold fsusie CV slice, or NULL when CV is disabled. diff --git a/man/GenotypeHandle-class.Rd b/man/GenotypeHandle-class.Rd index ce7fe80d..b7d25301 100644 --- a/man/GenotypeHandle-class.Rd +++ b/man/GenotypeHandle-class.Rd @@ -3,11 +3,108 @@ \docType{class} \name{GenotypeHandle-class} \alias{GenotypeHandle-class} +\alias{GenotypeHandle} \title{Genotype File Handle} +\usage{ +GenotypeHandle( + path = NULL, + plink1Prefix = NULL, + plink2Prefix = NULL, + bed = NULL, + bim = NULL, + fam = NULL, + pgen = NULL, + pvar = NULL, + psam = NULL, + ldMeta = NULL, + region = NULL, + genoMeta = NULL, + chroms = NULL, + format = NULL, + vcfArgs = list() +) +} +\arguments{ +\item{path}{Single file path (.vcf/.vcf.gz/.vcf.bgz/.bcf/.gds), or +\code{NULL}.} + +\item{plink1Prefix}{Path prefix for a PLINK1 triplet, or \code{NULL}.} + +\item{plink2Prefix}{Path prefix for a PLINK2 triplet, or \code{NULL}.} + +\item{bed, bim, fam}{Explicit paths to the PLINK1 triplet, or all \code{NULL}.} + +\item{pgen, pvar, psam}{Explicit paths to the PLINK2 triplet, or all +\code{NULL}.} + +\item{ldMeta}{Path to an LD-meta TSV file (columns \code{chrom}, +\code{start}, \code{end}, \code{path}; the \code{path} column may be +comma-separated as \code{ld_file,bim_file}). Requires \code{region}. The +constructor resolves the row covering \code{region}, then delegates to the +appropriate file-based handler. When the resolved row points at PLINK1 / +PLINK2 / VCF / GDS files, the corresponding reader is used; \code{.cor.xz} +(pre-computed LD-matrix) rows are not supported here -- use +\code{\link{loadLdMatrix}} for that case.} + +\item{region}{Region specification for \code{ldMeta} lookup: +\code{"chr:start-end"} string or a one-row data.frame with \code{chrom}, +\code{start}, \code{end}.} + +\item{genoMeta}{One-file-per-chromosome specification: a path to a +\code{#chr,path} meta file or a named character vector (names = +chromosomes, values = payload paths/prefixes). Optionally pass +\code{format} via \code{...} to force a single backend for every shard.} + +\item{chroms}{Optional character vector of chromosomes. With \code{genoMeta} +(a sharded, one-file-per-chromosome panel) only the shards whose chromosome +is listed are read, skipping the rest -- an I/O optimisation when the panel +is genome-wide but only a few chromosomes are needed. Chromosome labels are +compared canonically (\code{"chr1"}/\code{"1"} match, \code{23}/\code{X} +etc.). If none of the requested chromosomes are present in the panel, every +shard is read (so the caller's own absence check can report the mismatch). +Only meaningful with \code{genoMeta}; supplying it with any other source is +an error (a single-file panel has no per-chromosome shards to skip).} + +\item{vcfArgs}{Optional named list of arguments forwarded to +\code{VariantAnnotation::readVcf} when the source is a VCF.} +} +\value{ +A \code{GenotypeHandle} object. +} \description{ S4 container holding a path + format + metadata for lazy genotype access. Supports PLINK1 (.bed/.bim/.fam), PLINK2 (.pgen/.pvar/.psam), VCF/BCF, and GDS. + +Construct a \code{GenotypeHandle} from one of several input + forms. Exactly one of the following must be specified: + \describe{ + \item{\code{path}}{A single file path with a recognized extension: + \code{.vcf}, \code{.vcf.gz}, \code{.vcf.bgz}, \code{.bcf}, or + \code{.gds}. Format is auto-detected from the extension.} + \item{\code{plink1Prefix}}{A path prefix; the constructor appends + \code{.bed}, \code{.bim}, and \code{.fam} to locate the triplet.} + \item{\code{plink2Prefix}}{A path prefix; the constructor appends + \code{.pgen}, \code{.pvar} (or \code{.pvar.zst}), and \code{.psam} + to locate the triplet.} + \item{\code{bed} + \code{bim} + \code{fam}}{Explicit PLINK1 triplet. + The three files must share a stem; if they don't, use + \code{plink1Prefix} or arrange symlinks at a common stem.} + \item{\code{pgen} + \code{pvar} + \code{psam}}{Explicit PLINK2 + triplet. Same shared-stem requirement.} + \item{\code{genoMeta}}{One genotype file per chromosome. Either a path + to a whitespace/TSV meta file whose first column is the chromosome + (\code{#chr}) and second column the per-chromosome payload + (a \code{.bed}/\code{.pgen}/\code{.vcf[.gz]}/\code{.bcf}/\code{.gds} + file or a PLINK prefix; relative paths resolve against the meta + file's directory), or a named character vector mapping chromosome to + payload. All shards must share one format and identical sample IDs in + the same order. Extraction is routed to the correct file by the + requested region's chromosome.} + } + The constructor opens the file for metadata only (sample IDs and SNP + info); dosage extraction is deferred until \code{extractBlockGenotypes()} + is called. } \details{ Internal: the handle is the seed layer behind a genotype diff --git a/man/GenotypeHandle.Rd b/man/GenotypeHandle.Rd deleted file mode 100644 index 4562e8be..00000000 --- a/man/GenotypeHandle.Rd +++ /dev/null @@ -1,103 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/GenotypeHandle.R -\name{GenotypeHandle} -\alias{GenotypeHandle} -\title{Create a GenotypeHandle Object} -\usage{ -GenotypeHandle( - path = NULL, - plink1Prefix = NULL, - plink2Prefix = NULL, - bed = NULL, - bim = NULL, - fam = NULL, - pgen = NULL, - pvar = NULL, - psam = NULL, - ldMeta = NULL, - region = NULL, - genoMeta = NULL, - chroms = NULL, - format = NULL, - vcfArgs = list() -) -} -\arguments{ -\item{path}{Single file path (.vcf/.vcf.gz/.vcf.bgz/.bcf/.gds), or -\code{NULL}.} - -\item{plink1Prefix}{Path prefix for a PLINK1 triplet, or \code{NULL}.} - -\item{plink2Prefix}{Path prefix for a PLINK2 triplet, or \code{NULL}.} - -\item{bed, bim, fam}{Explicit paths to the PLINK1 triplet, or all \code{NULL}.} - -\item{pgen, pvar, psam}{Explicit paths to the PLINK2 triplet, or all -\code{NULL}.} - -\item{ldMeta}{Path to an LD-meta TSV file (columns \code{chrom}, -\code{start}, \code{end}, \code{path}; the \code{path} column may be -comma-separated as \code{ld_file,bim_file}). Requires \code{region}. The -constructor resolves the row covering \code{region}, then delegates to the -appropriate file-based handler. When the resolved row points at PLINK1 / -PLINK2 / VCF / GDS files, the corresponding reader is used; \code{.cor.xz} -(pre-computed LD-matrix) rows are not supported here -- use -\code{\link{loadLdMatrix}} for that case.} - -\item{region}{Region specification for \code{ldMeta} lookup: -\code{"chr:start-end"} string or a one-row data.frame with \code{chrom}, -\code{start}, \code{end}.} - -\item{genoMeta}{One-file-per-chromosome specification: a path to a -\code{#chr,path} meta file or a named character vector (names = -chromosomes, values = payload paths/prefixes). Optionally pass -\code{format} via \code{...} to force a single backend for every shard.} - -\item{chroms}{Optional character vector of chromosomes. With \code{genoMeta} -(a sharded, one-file-per-chromosome panel) only the shards whose chromosome -is listed are read, skipping the rest -- an I/O optimisation when the panel -is genome-wide but only a few chromosomes are needed. Chromosome labels are -compared canonically (\code{"chr1"}/\code{"1"} match, \code{23}/\code{X} -etc.). If none of the requested chromosomes are present in the panel, every -shard is read (so the caller's own absence check can report the mismatch). -Only meaningful with \code{genoMeta}; supplying it with any other source is -an error (a single-file panel has no per-chromosome shards to skip).} - -\item{vcfArgs}{Optional named list of arguments forwarded to -\code{VariantAnnotation::readVcf} when the source is a VCF.} -} -\value{ -A \code{GenotypeHandle} object. -} -\description{ -Construct a \code{GenotypeHandle} from one of several input - forms. Exactly one of the following must be specified: - \describe{ - \item{\code{path}}{A single file path with a recognized extension: - \code{.vcf}, \code{.vcf.gz}, \code{.vcf.bgz}, \code{.bcf}, or - \code{.gds}. Format is auto-detected from the extension.} - \item{\code{plink1Prefix}}{A path prefix; the constructor appends - \code{.bed}, \code{.bim}, and \code{.fam} to locate the triplet.} - \item{\code{plink2Prefix}}{A path prefix; the constructor appends - \code{.pgen}, \code{.pvar} (or \code{.pvar.zst}), and \code{.psam} - to locate the triplet.} - \item{\code{bed} + \code{bim} + \code{fam}}{Explicit PLINK1 triplet. - The three files must share a stem; if they don't, use - \code{plink1Prefix} or arrange symlinks at a common stem.} - \item{\code{pgen} + \code{pvar} + \code{psam}}{Explicit PLINK2 - triplet. Same shared-stem requirement.} - \item{\code{genoMeta}}{One genotype file per chromosome. Either a path - to a whitespace/TSV meta file whose first column is the chromosome - (\code{#chr}) and second column the per-chromosome payload - (a \code{.bed}/\code{.pgen}/\code{.vcf[.gz]}/\code{.bcf}/\code{.gds} - file or a PLINK prefix; relative paths resolve against the meta - file's directory), or a named character vector mapping chromosome to - payload. All shards must share one format and identical sample IDs in - the same order. Extraction is routed to the correct file by the - requested region's chromosome.} - } - The constructor opens the file for metadata only (sample IDs and SNP - info); dosage extraction is deferred until \code{extractBlockGenotypes()} - is called. -} -\keyword{internal} diff --git a/man/fsusieAffectedRegions.Rd b/man/fsusieAffectedRegions.Rd index 9a3ebb28..a46ab2ba 100644 --- a/man/fsusieAffectedRegions.Rd +++ b/man/fsusieAffectedRegions.Rd @@ -23,8 +23,8 @@ The sub-intervals of the functional grid where a credible set's band excludes zero (the genomic footprint of each fSuSiE effect), from the upstream \code{fsusieR::affected_reg}. Returns a \code{GRanges} with the credible-set label, its purity, and the effect \code{direction} -(\code{"pos"}/\code{"neg"}; upstream discards it). Requires an UNtrimmed fit; -empty otherwise. +(\code{"pos"}/\code{"neg"}; upstream discards it). Uses retained functional +summaries or an untrimmed fit; empty when those fields are absent. } \examples{ data(fsusieFineMappingExample) diff --git a/man/fsusieCredibleBand.Rd b/man/fsusieCredibleBand.Rd index 20a698bf..b1b371c8 100644 --- a/man/fsusieCredibleBand.Rd +++ b/man/fsusieCredibleBand.Rd @@ -23,8 +23,9 @@ For a functional-SuSiE (\code{fsusieR::susiF}) fine-mapping entry, returns the fitted effect curve and its credible band over the functional grid, one row per (credible set, grid position). Wraps the upstream fSuSiE band computation (fsusieR's internal \code{update_cal_credible_band.susiF} + -\code{get_fitted_effect}). Requires an UNtrimmed fit (the wavelet slots are -dropped by trimming); degrades to zero rows for non-fSuSiE or trimmed fits. +\code{get_fitted_effect}). Uses saved bands in trimmed fits, or computes +them from an untrimmed fit. Returns zero rows when the functional fields +are absent. } \examples{ data(fsusieFineMappingExample) diff --git a/pixi.toml b/pixi.toml index 3876de8f..a08cb005 100644 --- a/pixi.toml +++ b/pixi.toml @@ -9,7 +9,7 @@ platforms = [ [tasks] devtools_document = "cd $GITHUB_WORKSPACE; R -e 'devtools::document()'" -devtools_test = "cd $GITHUB_WORKSPACE; R -e 'devtools::test()'" +devtools_test = "cd $GITHUB_WORKSPACE; R -e 'devtools::test(stop_on_failure = TRUE)'" codecov = "cd $GITHUB_WORKSPACE; R -e 'covr::codecov(quiet = FALSE)'" rcmdcheck = "cd $GITHUB_WORKSPACE; R -e 'rcmdcheck::rcmdcheck()'" bioccheck_git_clone = "cd $GITHUB_WORKSPACE; R -e 'BiocCheck::BiocCheckGitClone()'" diff --git a/tests/testthat/test_fineMappingPipeline.R b/tests/testthat/test_fineMappingPipeline.R index d040a0b0..fb2dd94f 100644 --- a/tests/testthat/test_fineMappingPipeline.R +++ b/tests/testthat/test_fineMappingPipeline.R @@ -1923,10 +1923,16 @@ test_that("fineMappingPipeline(QtlDataset): jointSpec='trait' with fsusie wires ) ) expect_s4_class(res, "QtlFineMappingResult") - # fsusie over 2 traits emits one per-trait row each (shared functional fit). - expect_equal(nrow(res), 2L) + # fsusie fits ONE functional model across the traits and stores it once. + # Its curve spans the trait axis -- each trait is a POSITION on a shared + # grid -- so unlike mvsusie, whose posteriors carry a conditions axis + # that .fmSliceCond() slices per condition, there is nothing to + # partition: one row per trait would hold N copies of the same curve. + # The full membership stays visible in the `jointTraits` column. + expect_equal(nrow(res), 1L) expect_true(all(as.character(res$method) == "fsusie")) - expect_setequal(as.character(res$trait), c("ENSG_A", "ENSG_B")) + expect_equal(as.character(res$trait), "ENSG_A") + expect_equal(as.character(res$jointTraits), "ENSG_A;ENSG_B") expect_true("jointTraits" %in% pecotmr:::.tupleColumnNames(res)) }) @@ -1948,8 +1954,15 @@ test_that("fineMappingPipeline(QtlDataset): fsusie multi-trait per context dispa fineMappingPipeline(qd, methods = "fsusie", cisWindow = 1000L) ) expect_s4_class(res, "QtlFineMappingResult") - expect_equal(nrow(res), 2L) - expect_setequal(traitNames(res), c("ENSG_A", "ENSG_B")) + # fsusie fits ONE functional model across the traits and stores it once. + # Its curve spans the trait axis -- each trait is a POSITION on a shared + # grid -- so unlike mvsusie, whose posteriors carry a conditions axis + # that .fmSliceCond() slices per condition, there is nothing to + # partition: one row per trait would hold N copies of the same curve. + # The full membership stays visible in the `jointTraits` column. + expect_equal(nrow(res), 1L) + expect_equal(traitNames(res), "ENSG_A") + expect_equal(as.character(res$jointTraits), "ENSG_A;ENSG_B") expect_setequal(methodNames(res), "fsusie") }) @@ -3581,7 +3594,9 @@ test_that("fineMappingPipeline(QtlDataset): fsusie jointRegions=FALSE merges per region = regions, jointRegions = FALSE )) - expect_equal(nrow(res), 2L) + # One row for the joint fsusie fit (see above); the per-REGION merge is + # what `names(fit)` checks, and it is unaffected. + expect_equal(nrow(res), 1L) fit <- susieFit( res, studyName = "study1", diff --git a/tests/testthat/test_fineMappingRow.R b/tests/testthat/test_fineMappingRow.R index dea98f4c..dc5b1a03 100644 --- a/tests/testthat/test_fineMappingRow.R +++ b/tests/testthat/test_fineMappingRow.R @@ -972,7 +972,7 @@ test_that("credibleSetSummary aggregates across a collection with entry identity # test_fsusieAccessors.R when R/fsusieAccessors.R was folded in) # =========================================================================== -.fsa_makeFit <- function() { +.fsa_makeFit <- function(post_processing = "none") { set.seed(1) n <- 150L p <- 24L @@ -995,7 +995,7 @@ test_that("credibleSetSummary aggregates across a collection with entry identity Y = Y, pos = seq_len(J), L = 5, - post_processing = "none", + post_processing = post_processing, verbose = FALSE )) } @@ -1038,6 +1038,51 @@ test_that("fsusieCredibleBand returns a long effect + band table (lower <= effec expect_true(all(grepl("^fsusie_", cb$cs))) }) +test_that("trimmed fSuSiE keeps curves, bands and coloc inputs", { + skip_if_not_installed("fsusieR") + skip_if_not_installed("wavethresh") + for (mode in c("none", "TI")) { + fit <- .fsa_makeFit(mode) + expectedFit <- pecotmr:::.fsusiePopulateCredibleBand(fit) + tables <- list(list(sets = list(cs = list()))) + trimmed <- pecotmr:::trimFinemappingFit( + fit, seq_along(fit$alpha), "fsusie", tables + ) + entry <- unserialize(serialize(.fsa_entry(trimmed), NULL)) + saved <- entry@susieFit + expect_null(saved$fitted_wc) + expect_null(saved$fitted_wc2) + expect_identical(saved$fitted_func, fit$fitted_func) + expect_identical(saved$cred_band, expectedFit$cred_band) + expect_equal( + .fmrRowFsusieCredibleBand(entry), + .fmrRowFsusieCredibleBand(.fsa_entry(fit)) + ) + expect_equal( + .fmrRowFsusieAffectedRegions(entry), + .fmrRowFsusieAffectedRegions(.fsa_entry(fit)) + ) + # The five positional filter arguments were folded into + # ColocLbfFilterParam(); these are its defaults, pinned explicitly so + # a later default change cannot quietly alter what this asserts. + lbf <- pecotmr:::.colocExtractLbfFromEntry( + entry, + ColocLbfFilterParam( + filterLbfCs = FALSE, + secondary = NULL, + concentration = 0.5, + priorTol = 1e-9 + ) + ) + expected <- pecotmr:::.asLbfMatrix(fit) + colnames(expected) <- variantIds(entry) + expect_equal(lbf$lbf, expected) + colocResult <- coloc::coloc.bf_bf(lbf$lbf, lbf$lbf) + expect_true("PP.H4.abf" %in% names(colocResult$summary)) + } +}) + + test_that("fsusieAffectedRegions returns a GRanges with cs / purity / direction", { skip_if_not_installed("fsusieR") skip_if_not_installed("wavethresh") diff --git a/tests/testthat/test_jointEngine.R b/tests/testthat/test_jointEngine.R index 71e139ae..554ebce6 100644 --- a/tests/testthat/test_jointEngine.R +++ b/tests/testthat/test_jointEngine.R @@ -895,7 +895,7 @@ test_that(".runJointCell: cross-trait twas -> per-trait weight vectors", { expect_false(is.matrix(weights(pecotmr:::.collectionEntry(res, 1L)))) }) -test_that("fitJointGroup(Individual, Fm): fsusie returns one entry per trait", { +test_that("fitJointGroup(Individual, Fm): fsusie keeps one entry", { set.seed(12) n <- 10L X <- matrix( @@ -938,9 +938,15 @@ test_that("fitJointGroup(Individual, Fm): fsusie returns one entry per trait", { entries <- pecotmr:::fitJointGroup(grp, pipe, "fsusie", list()) expect_type(entries, "list") - expect_length(entries, 2L) # one per trait + expect_length(entries, 1L) expect_s4_class(entries[[1L]], "FineMappingRow") + expect_identical(captured$Y, Y) expect_equal(captured$pos, c(100, 200)) # functional domain threaded + # The trait list lives in the collection's `jointTraits` column, so the + # fit carries only the grid POSITIONS -- the one piece that column + # cannot express. + expect_null(entries[[1L]]@susieFit$trait_names) + expect_equal(entries[[1L]]@susieFit$trait_positions, c(100, 200)) }) test_that("fitJointGroup(Individual, Fm): fsusie without pos errors; unknown token errors", { @@ -1900,7 +1906,7 @@ test_that("fitJointGroup(Individual, Fm): fsusie honest per-fold CV is attached" ) entries <- pecotmr:::fitJointGroup(grp, pipe, "fsusie", list()) expect_true(cvCalled) # CV path exercised - expect_length(entries, 2L) + expect_length(entries, 1L) }) test_that("fitJointGroup(Individual, Fm): SER pre-screen skips when < 2 survivors", {