NEWS.md
The optimal weights are computed with an accurate quadrature and a longer Hermite expansion. The moments behind Optimal_Weights_M() (the _N and _sparse schemes of GLOW_Fisher and GLOW_SKAT) are one-dimensional integrals of the transformation against the normal density. Up to 0.3.3 they were computed by integrate() at its default relative tolerance (), and Mehler’s expansion of the covariances stopped at order 8. On the alternative scale the Fisher transform is stretched to with up to 8 or more, its order-2 coefficient grows to about 100, and the adaptive rule returned it with an error of while reporting ; the covariance of a window with near-duplicate variants (condition number near ) turned that into weights off by 9 percent of the largest weight. Reported by Victoria Pan (WPI, 2026-10-29) with a reproducible script. Now every integral uses a composite 10-point Gauss-Legendre rule on with panels of width 0.1, cut where the argument of is zero (the kink of the two-sided transforms), the Hermite functions come from their recurrence to order 64 (ORD = 1:64; hermite() accepts any positive order), and each distinct (mean, scale) pair is integrated once per call. Against a direct integration of every pairwise covariance on a 20-variant window with a near-duplicate pair, the weights of the normal approximation now agree to (APR_N) and (BE_N) of the largest weight, the sparse weights to (APR_sparse) and (BE_sparse); the order-8 series with accurate coefficients was within and for the two _N schemes (tests/testthat/test-optimal-weights-quadrature.R, with the reference regenerated by tests/slow/optimal-weights-direct-reference.R). The sparse covariance keeps the exact mixture variance on its diagonal, now by the law of total variance with centered variances, and the normal covariance keeps the truncated variance, which measured better for its weights. The null covariance cached by CovMT_mix() is keyed by its numerical inputs (the null coefficients, M and the orders), so two closures with the same body and different captured values no longer share an entry. The function is also faster: 0.15 s against 0.39 s for 20 variants, 0.71 s against 2.19 s for 100 and 10.3 s against 25.0 s for 1,000, on one core. What changes for users: the _N and _sparse weights of GLOW_Fisher, and through them its optimal-weight p-values and the omnibus columns; GLOW_SKAT’s weights change only by numerical error (its expansion is exact at order two); the equal-weight tests, Burden and the single-variant scan are unchanged. The quadrature is validated for the Fisher and polynomial transforms; a custom g with a kink away from zero, a cusp, a discontinuity or rapid growth needs its own convergence check (see ?Optimal_Weights_M).
glow_region_test(debug_errors = TRUE) names the failing call and the stack. When a region’s test fails, the message of its status = "error" result carries, after the error’s own message, the call that signaled it and the names of the functions from glow_region_test() down to the signal, so the failure can be located without rerunning the region. Off by default; results unchanged. Suggested by Alberto Brusati.
Optimal_Weights_M() (the wts_BE_N and wts_APE_N schemes) is computed on the standardized scale. Under the alternative the Z-scores are approximated by , whose covariance has a diagonal above one. The moment routines are written for a standard scale: get_r_tilde() takes the standard deviation of each , and CovM_gXgY() applies Mehler’s expansion, an expansion in a correlation whose entries are bounded by one, which is what lets the sum stop at order 8. Up to 0.3.2, as in the legacy code, Optimal_Weights_M() passed the covariance M + MVM into the expansion and its diagonal, a variance, as the standard deviation, so the order-8 term grew with the entries of M + MVM instead of shrinking (on the ALS chromosome 21 windows of the report, up to times the bound). Reported by Victoria Pan (WPI, 2026-10-23) with a reproducible script. Now Optimal_Weights_M() writes with and , get_r_tilde() receives , and CovM_gXgY() receives the correlation matrix with the scales through its new SD1 and SD2 arguments (the Hermite coefficients are those of ); it refuses a matrix without a unit diagonal. Its identity fast path now uses both mean vectors. For the result matches the closed forms (tests/testthat/test-optimal-weights-normal-approx.R). The _N weights change wherever MVM is not zero, that is wherever some PI lies strictly between 0 and 1; the _sparse and Burden weights and every test statistic are unchanged. The p-values of the optimal-weight tests also change through GFisher 0.3.2, whose HYB method applied unequal weights twice (its NEWS); the equal-weight tests GLOW_SKAT and GLOW_Fisher are unaffected by either.GenomeInfoDb. Since Bioconductor 3.22 the seqnames() generic lives in the new Seqinfo package, and GenomicFeatures and GenomicRanges no longer install GenomeInfoDb. On a current R the two promoter masks therefore stopped with “The promoter masks need the GenomeInfoDb package” although everything they use was installed, which is how the public check of GLOWpipeline 0.2.1 failed. The one call now goes through GenomicRanges::seqnames(), which both Bioconductor generations export, and GenomeInfoDb leaves Suggests. The masks need GenomicFeatures and TxDb.Hsapiens.UCSC.hg38.knownGene (with their own dependencies) only. No result changes.citation("GLOWr") and the README) now cites Zhang, Liu, Landers and Wu (2026), The Annals of Applied Statistics 20(3), 2074-2097, doi:10.1214/26-AOAS2223, in place of “in revision”. Documentation only; no code change.define_regions_noncoding() builds, once per chromosome and mask, the regions of STAARpipeline’s eight gene-centric non-coding masks (upstream, downstream, UTR, promoter_CAGE, promoter_DHS, enhancer_CAGE, enhancer_DHS, ncRNA). A region is one gene. Its candidates are the records that the annotation links to the gene wherever they lie on the chromosome, stored as a list column variant_rows of row numbers in the chromosome’s variant index. extract_variant_set() takes such a region in place of its interval lookup, after validating the list (unique, increasing, within the index). Everything after the lookup is unchanged. noncoding_filter(mask, filter_spec) returns the per-set filter (the run’s filter with no annotation clause), and noncoding_masks() the eight names. The masks read the aGDS fields GENCODE.Category, GENCODE.Info, GeneHancer, CAGE and DHS of the bundled catalog, which annotate_favor() writes by default, in both the FAVOR v1 and the FAVOR 2.0 layouts. The gene is read by the meaning of the field: details in parentheses are removed, ; separates the groups of a multi-category record, , separates genes, NONE is a sentinel, and hyphenated names stay whole. A record annotated upstream;downstream therefore joins both masks, a UTR record naming two genes joins both sets, and an ncRNA record with several groups links them all. STAARpipeline’s literal string rules are available as rule = "staar", and they gave the same candidate sets variant for variant on chromosome 22 of the 1000 Genomes Project. The two promoter masks keep the records within 3 kb of a transcription start site of TxDb.Hsapiens.UCSC.hg38.knownGene. That package is now suggested, together with GenomicFeatures, GenomicRanges, IRanges and GenomeInfoDb.load_B_model() checks the file’s format, not the package version that saved it. The loader warned whenever the model had been saved by a different GLOWr version, so a valid model file produced a warning on every load after the next release. It now compares the envelope’s format_version with the format this version reads (1), warns when the stamp differs, is absent or is not a single number, and keeps the saving version in the file’s metadata as provenance. A model saved by an earlier GLOWr release carries format 1 in its envelope and loads without a warning and predicts as before (the model sent to Alberto Brusati on 2026-09-01, saved by 0.1.0, is one).Changes for whole-genome-sequencing scans, after Alberto Brusati’s Terra runs of chromosomes with eleven million variant records. Results are unchanged wherever the previous version completed.
build_variant_index() reads a GDS file’s chromosome, position and id vectors once. extract_variant_set() and glow_region_test() take it as variant_index and find a region’s records by binary search, setting the GDS filter by file position, instead of re-reading the chromosome’s variant table and matching ids for every region. On a chromosome of one million records a 2 kb window cost 0.40 s in the lookup before and about 0.02 s after; the saving grows with the chromosome’s record count. A call without an index builds one for itself and selects exactly as before. count_index_records() counts the records of an index inside genomic spans.filter_variants_ld() computes the correlation of its NA-free genotype matrix as one BLAS cross-product of the standardized columns instead of cor(use = "pairwise.complete.obs"), and its greedy prune keeps each variant’s count of above-threshold partners up to date instead of recomputing the whole threshold mask after every removal. The selection rule and the tie-breaking are unchanged, so the kept variants are the same. A gene of 2,018 rare SNVs took 169 s in this step before. A matrix with missing values keeps the previous path.filter_variants_ld() drops the columns that are exact linear combinations of the columns before them. It did this with LINPACK qr(), which cycles every dropped column to the end of the matrix by shifting the remaining columns, so a gene with far more rare variants than samples (14,038 columns, 3,195 kept, at 3,202 samples) spent 367 s there. The step is now a blocked in-order Gram-Schmidt with qr()’s rule and tolerance (a column is dropped when its residual against the kept columns before it is below 1e-07 of its norm), in BLAS calls, stopping once the kept count reaches the number of samples carrying a variant. The kept columns are the same, checked on every chromosome 22 gene of the 1000 Genomes cohort; that gene’s step takes a few seconds. A matrix with missing values now stops with a message, where qr() stopped with an opaque one.SPAtest::ScoreTest_SPA() skips a column whose min(sum(g), sum(2 - g)) is below its minimum, a guard written for dosages in [0, 2]. A collapsed burden column (the row sum of many ultra-rare variants) with a mean above 1 fails it, came back NA, and glow_test() stopped on the gene. SPAtest’s own computation uses the normal approximation for such a column past the guard, so getZ_marg_score_binary_SPA() now gives it the standard score Z, warns with the count, and returns the count as n_spa_fallback. This matches the GLOW methodology paper, which applied SPA to the non-collapsed variants only.marginal_scan() gains n_cores (default 1): the variant chunks are spread over forked workers, each with its own GDS handle, and the result table is the same as with one core. Not available on Windows.vignette("favor-annotation-quick-start") walks a new user from when to use annotate_favor(), through what to have ready and which settings to choose, to what the call returns and how to check it. It points into the full vignette for the evidence, the cost table and the matching procedure, which are unchanged. Written from user feedback. No code changed.SeqArray::seqVCF2GDS() stores a VCF FILTER of "." (filters not applied) as NA in annotation/filter. The variant filter compared that value with "PASS", obtained NA, and the region stopped with “missing value where TRUE/FALSE needed”. A missing QC value now counts as not passing, so such variants are excluded like any other non-passing variant, and ?variant_filter says so. A GDS converted from PLINK carries PASS for every variant and is unaffected. Reported by Alberto Brusati on an aGDS built from a VCF.evaluate_PI_models() and plot_PI_roc() score tied predictions correctly: the AUC is the rank-based statistic with ties at half credit, and the ROC curve takes one point per distinct threshold. Before, an ensemble member whose LASSO kept no feature (a constant prediction) scored an AUC of 1.0 whenever the cases were listed first, which is what the earlier ALS evaluations’ “AUC = 1” models were.annotate_favor() is rewritten around a corrected matching rule and gains a second database backend. The returned table, the aGDS nodes and the provenance change shape, hence the minor version. The arguments of 0.1.x keep their names and positions.
match_method = "flexible" the lookup key of an SNV is normalized against FAVOR’s reference base, by swapping REF and ALT or complementing both alleles. Every input row records its outcome in match_tier (exact, swapped, flipped, flipped_swapped, position_only, or a reason: rsid_conflict, unmatched_alt, unmatched_ref, uncovered, no_database, unsupported_allele) and the database row it used in favor_key. The 0.1.x flexible matching had no strand tier: on genotyping-chip data half the SNVs fell through to the mean over the three alternative alleles at their position and stored a value that belongs to no variant. A matched variant now always carries its own row’s values. Indels match by the exact key only. The cohort’s alleles, genotypes and keys are never modified.annotation/id, or an rsID column of a data.frame input) is compared with the FAVOR row’s rsID and recorded in rsid_check (same, differs, chip_none, favor_none). The new argument rsid_policy is "require" by default: a swapped or flipped match that the rsID does not confirm is withheld as rsid_conflict, so a rewritten key is accepted only when dbSNP confirms it or the input carries no rsID to check. Exact matches are never withheld. "record" keeps every match and only records the check. The rsID is never used to find a row.aPC.Protein now names apc_protein_function_v3, aPC.Conservation names apc_conservation_v2 and aPC.LocalDiversity names apc_local_nucleotide_diversity_v3, the versioned nodes the FAVOR databases carry.favor_db_format = "parquet" (detected under the default "auto") reads the FAVOR team’s per-chromosome Parquet files directly with the arrow package (now in Suggests), one row group at a time and only the leaf columns requested. It serializes GeneHancer and the GENCODE information fields into the v1 string layout. Both backends return identical outcomes for identical content. A field the database lacks is skipped with a message. FAVOR 2.0 has no apc_conservation, apc_local_nucleotide_diversity, apc_local_nucleotide_diversity_v2 or apc_proximity_to_coding.CHR-POS-NA-NA) are averaged over the SNV rows at the position only; a position with no SNV row is uncovered.annotation/info/favor_match_tier and annotation/info/favor_rsid_check beside the FunctionalAnnotation folder, whose attributes record the database format and files (with their SHA-256 when a SHA256SUMS file sits beside them), the favor_release label, the GLOWr version, the matching settings, the rsID policy and source, and the counts per outcome and per check. The writers align rows to the GDS through variant.id.favor_db_format, favor_release and rsid_policy, placed after verbose so that positional calls written for 0.1.x keep their meaning. GDS input returns variant_id and rsID beside VarInfo. A VarInfo position must be a positive integer in decimal digits. A fractional or zero position used to be truncated onto another coordinate. The CSV backend fetches the rows at the input positions by an xsv join on position, or by fread() with a filter. xsv file paths are shell-quoted.favor-annotation explains the matching rule, what a tier claims, the rsID evidence, position-only keys, the optional PLINK normalization of chip data and the measured cost. inst/scripts/annotate_favor_batch.R exposes --favor-db-format, --favor-release and --rsid-policy.annotate_favor()’s xsv fast path (exact matching): variants falling in the second and later FAVOR chunk files of a chromosome silently lost their annotations (the per-chunk left-joins made the downstream deduplication keep an empty row instead of the real match). The per-chunk joins are now inner joins, and every xsv invocation’s exit status is checked, so a failed join aborts instead of writing a silently under-annotated aGDS.annotate_favor()’s default features now covers the full FAVOR Essential DB annotation content (30 columns, previously the 20 purely numeric ones), so an aGDS built with the defaults carries the categorical features (e.g. GENCODE categories) that the built-in coding masks filter on. String features are first-class on the R join path: flexible-match aggregation is type-aware (numeric features average across matches; character features take the first non-empty value), and na_handling = "zero" no longer writes 0 into character columns.annotation_predicate() filter grammar for variant_filter(): a clause condition may now test nonempty/empty, in/not_in, or numeric gt/ge/lt/le alongside the existing set-membership conditions. Missing annotation values never satisfy a positive test.calibrate_pvalues() gains inflation_only: when TRUE, a calibration factor below 1 is clamped to 1, reproducing the field-standard one-sided genomic-control convention (correct inflation, leave deflation alone). Default FALSE keeps the previous two-sided behavior.src/Makevars and src/Makevars.win); the Rtools toolchain does not auto-link them, so this fixes the Windows R CMD check.PI is described throughout as the variant-importance score, and the methodology reference now cites the GLOW methods paper (Zhang, Liu, Landers, and Wu, Annals of Applied Statistics, in revision).Initial public release.
B and the variant-importance score PI, and their combination into per-variant weights used by the set-level tests.