add and clean up vignettes

NAMESPACE

@@ -42,6 +42,8 @@ export(bayesCWeights)
 export(bayesLWeights)
 export(bayesNWeights)
 export(bayesRWeights)
+export(buildLdEigen)
+export(buildLdScore)
 export(buildMrmashPriorMatrices)
 export(buildTopLoci)
 export(calculateFeatureScores)
@@ -84,7 +86,6 @@ export(filterInvalidSummaryStat)
 export(filterMixtureComponents)
 export(filterRelatedness)
 export(filterVariantsByLdReference)
-export(findOverlappingRegions)
 export(fineMappingPipeline)
 export(fineMappingRow)
 export(finemapCtwasRegions)
@@ -177,14 +178,13 @@ export(getSnpIdx)
 export(getStandardized)
 export(getStudy)
 export(getSumStats)
-export(getSumstatDf)
+export(getSumStatsDf)
 export(getSusieAlpha)
 export(getSusieFit)
 export(getSusieResult)
 export(getTauBlocks)
 export(getTopLoci)
 export(getTraitName)
-export(getTraitNames)
 export(getTraitPosition)
 export(getTraitRun)
 export(getTraitRuns)
@@ -210,15 +210,12 @@ export(lassosumRss)
 export(lassosumRssWeights)
 export(lbfToAlpha)
 export(ldClumpByScore)
-export(ldLoader)
 export(ldMismatchQc)
 export(ldPruneByCorrelation)
 export(learnTwasWeights)
 export(loadGenotypeRegion)
 export(loadGwasSumStatsFromManifest)
-export(loadLdBlock)
 export(loadLdMatrix)
-export(loadLdSketch)
 export(loadMultiStudyQtlDatasetFromManifest)
 export(loadQtlDatasetFromManifest)
 export(loadQtlSumStatsFromManifest)
@@ -240,7 +237,7 @@ export(mergeSusieCs)
 export(mergeVariantInfo)
 export(metaAnalysisPerCondition)
 export(metaSldscRandom)
-export(mrAshRssWeights)
+export(mrashRssWeights)
 export(mrashWeights)
 export(mrmashRssWeights)
 export(mrmashWeights)
@@ -271,7 +268,6 @@ export(readSldscAnnot)
 export(readSldscFrq)
 export(readSldscTrait)
 export(regionToDf)
-export(regionsOverlap)
 export(resolveWeights)
 export(sanitizeMashData)
 export(scadRssWeights)
@@ -300,9 +296,8 @@ export(twasWeightsPipeline)
 export(twasWeightsRow)
 export(twasZ)
 export(updateMashModelCov)
-export(variantIdToDf)
 export(waldTestPval)
-export(writeSumstatsVcf)
+export(writeSumStatsVcf)
 exportClasses(AnnotationMatrix)
 exportClasses(ColocBoostResult)
 exportClasses(ColocResult)
@@ -416,13 +411,12 @@ exportMethods(getSnpIdx)
 exportMethods(getStandardized)
 exportMethods(getStudy)
 exportMethods(getSumStats)
-exportMethods(getSumstatDf)
+exportMethods(getSumStatsDf)
 exportMethods(getSusieAlpha)
 exportMethods(getSusieFit)
 exportMethods(getTauBlocks)
 exportMethods(getTopLoci)
 exportMethods(getTraitName)
-exportMethods(getTraitNames)
 exportMethods(getTraitPosition)
 exportMethods(getTraitRun)
 exportMethods(getTraitRuns)
@@ -449,7 +443,7 @@ exportMethods(show)
 exportMethods(subsetChr)
 exportMethods(subsetRegion)
 exportMethods(twasWeightsPipeline)
-exportMethods(writeSumstatsVcf)
+exportMethods(writeSumStatsVcf)
 importClassesFrom(GenomicRanges,GRanges)
 importClassesFrom(MultiAssayExperiment,MultiAssayExperiment)
 importClassesFrom(S4Vectors,DataFrame)

R/AllClasses.R

@@ -85,7 +85,7 @@ setClassUnion("LdMixtureWeights", c("numeric", "NULL"))
 #
 # getZ / getN / getMaf / nSnps are
 # defined once on SumStatsBase (they only delegate to getSumStats); subsetChr /
-# getVarY / getSumStats / getSumstatDf stay on the concrete subclass because
+# getVarY / getSumStats / getSumStatsDf stay on the concrete subclass because
 # they rely on the tuple shape (3-tuple QtlSumStats, 1-tuple GwasSumStats).
 # =============================================================================
 

R/AllGenerics.R

@@ -23,11 +23,30 @@ NULL
 # =============================================================================
 
 #' @title Estimate SNP Heritability
-#' @description Estimate SNP heritability from GWAS summary statistics using one
-#'   of three methods: LDER, g-LDSC, or HDL/sHDL.
+#' @description Estimate SNP heritability from GWAS summary statistics using
+#'   one of four methods: S-LDSC, g-LDSC, LDER, or HDL/sHDL. Supplying
+#'   \code{annotations} makes any of them stratified.
 #' @param sumstats A \code{GwasSumStats} object.
-#' @param ldRef An \code{LdStatistic} object (method-appropriate subclass).
-#' @param method Character, one of "lder", "gldsc", "hdl".
+#' @param ldRef An \code{LdStatistic} object (method-appropriate subclass):
+#'   an \code{\link{LdScore}} for \code{"sldsc"} / \code{"gldsc"}, an
+#'   \code{\link{LdEigen}} for \code{"lder"} / \code{"hdl"}. Build either
+#'   from loaded LD with \code{\link{buildLdScore}} /
+#'   \code{\link{buildLdEigen}}.
+#' @param method Character, one of \code{"lder"} (the default),
+#'   \code{"sldsc"}, \code{"gldsc"} or \code{"hdl"}:
+#'   \describe{
+#'     \item{\code{"sldsc"}}{Stratified LD score regression (Finucane et al.
+#'       2015): weighted least squares of \eqn{\chi^2} on LD scores. Needs
+#'       only per-variant scores when unstratified, so it is the cheapest of
+#'       the four and the only one that runs on a scores-only
+#'       \code{LdScore}.}
+#'     \item{\code{"gldsc"}}{Generalized LD score regression (Xiong et al.
+#'       2024): GLS using the full per-block LD matrices.}
+#'     \item{\code{"lder"}}{LD eigenvalue regression (Song et al. 2022).}
+#'     \item{\code{"hdl"}}{High-definition likelihood (Ning et al. 2020),
+#'       stratified as sHDL. Models reference-panel noise explicitly, so it
+#'       shrinks hard when \code{nRef} is small.}
+#'   }
 #' @param annotations An \code{AnnotationMatrix} object, or NULL for
 #'   unstratified estimation.
 #' @param local Logical, whether to compute per-block local estimates.
@@ -36,17 +55,23 @@ NULL
 #'   this study; \code{NULL} matches all studies.
 #' @return An \code{H2Estimate} object.
 #' @examples
-#' data(ldEigenExample)
-#' gr <- GenomicRanges::GRanges("chr1",
-#'   IRanges::IRanges(seq(50, by = 100, length.out = 20), width = 1))
-#' S4Vectors::mcols(gr) <- S4Vectors::DataFrame(SNP = paste0("rs", 1:20),
-#'   A1 = "A", A2 = "G", Z = rnorm(20), N = 10000L)
+#' data(ldScoreExample)
+#' # The estimators index the z-scores by the reference's own variant order,
+#' # so build the sumstats from the reference rather than alongside it.
+#' gr <- as(ldScoreExample, "GRanges")
+#' set.seed(1)
+#' S4Vectors::mcols(gr) <- S4Vectors::DataFrame(
+#'   SNP = names(ldScoreExample),
+#'   A1 = S4Vectors::mcols(ldScoreExample)$A1,
+#'   A2 = S4Vectors::mcols(ldScoreExample)$A2,
+#'   Z = rnorm(length(ldScoreExample)),
+#'   N = 100000L)
 #' panel <- readGenotypes(
 #'   system.file("extdata", "toy_ref.bed", package = "pecotmr")
 #' )
 #' ss <- GwasSumStats(study = "trait1", entry = list(gr),
-#'   genome = "hg19", ldSketch = panel)
-#' estimateH2(ss, ldEigenExample, method = "lder")
+#'   genome = getGenome(ldScoreExample), ldSketch = panel)
+#' estimateH2(ss, ldScoreExample, method = "sldsc")
 #' @export
 setGeneric(
     "estimateH2",
@@ -561,9 +586,9 @@ setGeneric("getSumStats", function(x, ...) standardGeneric("getSumStats"))
 #' @return A \code{data.frame}.
 #' @examples
 #' data(qtlSumStatsExample)
-#' getSumstatDf(qtlSumStatsExample)
+#' getSumStatsDf(qtlSumStatsExample)
 #' @export
-setGeneric("getSumstatDf", function(x, ...) standardGeneric("getSumstatDf"))
+setGeneric("getSumStatsDf", function(x, ...) standardGeneric("getSumStatsDf"))
 
 #' @title Get the Embedded QtlDataset List
 #' @description Return the named list of \code{QtlDataset} objects carried by a
@@ -1203,8 +1228,10 @@ setGeneric("getDataType", function(x, ...) standardGeneric("getDataType"))
 #' Write summary statistics or fine-mapping results to VCF/BCF
 #'
 #' Creates a VCF object from GWAS summary statistics or fine-mapping results and
-#' writes it to disk. Supports bgzipped VCF (.vcf.gz/.vcf.bgz) and BCF (.bcf)
-#' output formats via VariantAnnotation and Rsamtools.
+#' writes it to disk via VariantAnnotation and Rsamtools. Supports bgzipped
+#' VCF (.vcf.gz/.vcf.bgz); BCF (.bcf) additionally needs a working
+#' \code{Rsamtools::asBcf()}, which recent Rsamtools releases disable, so
+#' prefer .vcf.bgz.
 #'
 #' @param x Input data: a \code{GwasSumStats} object, a \code{FineMappingResult}
 #'   object, or a data.frame with columns \code{chrom}, \code{pos}, \code{ref},
@@ -1230,11 +1257,11 @@ setGeneric("getDataType", function(x, ...) standardGeneric("getDataType"))
 #' @return Invisible path to the written file.
 #' @examples
 #' data(gwasSumStatsS4Example)
-#' writeSumstatsVcf(
+#' writeSumStatsVcf(
 #'   gwasSumStatsS4Example, outputPath = tempfile(fileext = ".vcf"))
 #' @export
-setGeneric("writeSumstatsVcf", function(x, outputPath, sampleName = NULL, ...) {
-    standardGeneric("writeSumstatsVcf")
+setGeneric("writeSumStatsVcf", function(x, outputPath, sampleName = NULL, ...) {
+    standardGeneric("writeSumStatsVcf")
 })
 
 # =============================================================================
@@ -1262,7 +1289,8 @@ setGeneric("getContexts", function(x) standardGeneric("getContexts"))
 
 #' @title Get Unique Trait Names
 #' @description Return the unique trait identifiers carried by a collection
-#'   class (e.g., \code{QtlSumStats}).
+#'   class (e.g., \code{QtlSumStats}), or the trait names of an
+#'   \code{\link{SldscData}}.
 #' @param x The object.
 #' @return Character vector of unique trait names.
 #' @examples
@@ -1893,38 +1921,6 @@ setGeneric("getFrqData", function(x) standardGeneric("getFrqData"))
 #' @export
 setGeneric("getTraitRuns", function(x) standardGeneric("getTraitRuns"))
 
-#' @title Get the trait names from an SldscData
-#' @param x An \code{\link{SldscData}} object.
-#' @return A character vector of trait names.
-#' @rdname getTraitNames
-#' @examples
-#' mkRun <- function(cats) {
-#'   n <- length(cats)
-#'   list(categories = cats, tau = setNames(rep(1e-7, n), cats),
-#'     tauSe = setNames(rep(3e-8, n), cats),
-#'     enrichment = setNames(rep(2, n), cats),
-#'     enrichmentSe = setNames(rep(0.4, n), cats),
-#'     enrichmentP = setNames(rep(0.01, n), cats),
-#'     propH2 = setNames(rep(0.2, n), cats),
-#'     propSnps = setNames(rep(0.1, n), cats), h2g = 0.3,
-#'     tauBlocks = matrix(1e-7, 10, n, dimnames = list(NULL, cats)),
-#'     nBlocks = 10L)
-#' }
-#' annot <- data.frame(CHR = c(1, 1, 1, 2, 2, 2), SNP = paste0("rs", 1:6),
-#'   annot_A = c(1, 0, 1, 0, 1, 0), annot_B = c(2.1, 1.8, 2.5, 1.9, 2.3, 2))
-#' frq <- data.frame(CHR = c(1, 1, 1, 2, 2, 2), SNP = paste0("rs", 1:6),
-#'   MAF = rep(0.2, 6))
-#' mkTrait <- function() {
-#'   list(single = list(mkRun(c("annot_A_0", "baselineLD_0")),
-#'     mkRun(c("annot_B_0", "baselineLD_0"))),
-#'     joint = mkRun(c("annot_A_0", "annot_B_0", "baselineLD_0")))
-#' }
-#' traits <- setNames(list(mkTrait(), mkTrait()), c("traitX", "traitY"))
-#' sd <- SldscData(annot = annot, frq = frq, traits = traits)
-#' getTraitNames(sd)
-#' @export
-setGeneric("getTraitNames", function(x) standardGeneric("getTraitNames"))
-
 #' @title Get the annotation column names from an SldscData
 #' @param x An \code{\link{SldscData}} object.
 #' @return A character vector of annotation column names.

R/GenotypeHandle.R

@@ -310,7 +310,7 @@ GenotypeHandle <- function(
         msg <- glue(
             "GenotypeHandle: the LD-meta row for region {regionStr} points ",
             "at a pre-computed correlation matrix ({ldPath}). Use ",
-            "loadLdMatrix() / loadLdSketch() for .cor.xz inputs; ",
+            "loadLdMatrix() for .cor.xz inputs; ",
             "GenotypeHandle accepts only genotype payloads (VCF/GDS/PLINK)."
         )
         abort(msg)

R/LdEigen.R

@@ -110,6 +110,103 @@ LdEigen <- function(
     obj
 }
 
+#' @title Build an LdEigen from loaded LD
+#' @description Eigendecompose already-loaded LD, block by block, into the
+#'   \code{LdEigen} that \code{\link{estimateH2}} consumes for
+#'   \code{method = "lder"} and \code{method = "hdl"}. This is the supported
+#'   route from an LD reference to an h2 input: \code{\link{LdEigen}} itself
+#'   is the low-level constructor and expects the decompositions to have been
+#'   computed already.
+#'
+#'   No file I/O happens here. Read the LD first with
+#'   \code{\link{loadLdMatrix}} -- passing a vector of regions returns one
+#'   \code{LdData} per block -- then pass the result in.
+#' @param ldBlockData An \code{\link{LdData}}, or a list of them, covering
+#'   the reference variants. An \code{LdData} whose correlation is a list of
+#'   per-block matrices contributes one LD block per matrix; one whose
+#'   correlation is a single matrix contributes one block. Blocks are kept in
+#'   the order given, and that order defines the variant order of the result.
+#'
+#'   Block structure matters downstream: \code{\link{estimateH2}} takes its
+#'   standard error from a delete-one-block jackknife, so it requires at
+#'   least two blocks and wants many more. Loading a whole region as one
+#'   dense matrix yields a single block.
+#' @param nRef Integer, the LD reference panel sample size. Defaults to the
+#'   size recorded by the supplied \code{LdData}; required when they record
+#'   none, and an error when they disagree.
+#' @param inSample Logical, whether the reference is the GWAS cohort itself.
+#'   Selects LDER's weighting and HDL's finite-reference correction.
+#' @param genome Character, genome build. Defaults to the build recorded by
+#'   the supplied \code{LdData}, which \code{\link{loadLdMatrix}} leaves
+#'   unset.
+#' @param eigenvalueTruncation Numeric in (0, 1]; the proportion of each
+#'   block's eigenvalue mass to retain. \code{1} (the default) keeps every
+#'   component; HDL conventionally uses \code{0.9}.
+#' @return An \code{LdEigen} over every variant in \code{ldBlockData}.
+#' @seealso \code{\link{buildLdScore}} for the g-LDSC input,
+#'   \code{\link{estimateH2}}, \code{\link{loadLdMatrix}}
+#' @examples
+#' meta <- system.file("extdata", "ld_reference", "ld_meta_file.tsv",
+#'   package = "pecotmr")
+#' ld <- loadLdMatrix(meta, region = "chr22:10000000-19000000")
+#' ldEigen <- buildLdEigen(ld, genome = "hg38")
+#' ldEigen
+#' length(getEigenList(ldEigen))
+#' @export
+buildLdEigen <- function(
+    ldBlockData,
+    nRef = NULL,
+    inSample = FALSE,
+    genome = NA_character_,
+    eigenvalueTruncation = 1
+) {
+    prep <- .ldRefPrepare(ldBlockData, nRef, genome)
+    eigenList <- map2(
+        prep$blocks,
+        prep$snpIdx,
+        .ldEigenOneBlock,
+        truncation = eigenvalueTruncation
+    )
+    LdEigen(
+        snpInfo = prep$snpInfo,
+        eigenList = eigenList,
+        ldBlocks = prep$ldBlocks,
+        nRef = prep$nRef,
+        inSample = inSample,
+        genome = prep$genome,
+        eigenvalueTruncation = eigenvalueTruncation
+    )
+}
+
+# One block's decomposition, truncated to the leading components carrying
+# `truncation` of its eigenvalue mass.
+# @noRd
+.ldEigenOneBlock <- function(block, snpIdx, truncation) {
+    e <- eigen(block$R, symmetric = TRUE)
+    keep <- .ldEigenKeep(e$values, truncation)
+    list(
+        values = e$values[keep],
+        vectors = e$vectors[, keep, drop = FALSE],
+        snpIdx = as.integer(snpIdx)
+    )
+}
+
+# Rounding can leave a near-singular block with slightly negative trailing
+# eigenvalues, so the cumulative mass is taken over the clamped values --
+# otherwise the running total can dip and pick the wrong cut point.
+# @noRd
+.ldEigenKeep <- function(values, truncation) {
+    if (truncation >= 1) {
+        return(seq_along(values))
+    }
+    positive <- pmax(values, 0)
+    total <- sum(positive)
+    if (total <= 0) {
+        return(seq_along(values))
+    }
+    seq_len(which(cumsum(positive) / total >= truncation)[[1L]])
+}
+
 #' @describeIn LdEigen-class Refused. Subsetting would narrow the variants
 #'   while \code{eigenList} -- a slot, because it is per block rather than
 #'   per variant -- stayed as it was, leaving decompositions that describe