Repository navigation
Conversation
TWAS weights released on Synapse use the legacy list layout, which pecotmr 0.8.x does not read. This adds a protocol-level conversion layer, as agreed, rather than a backwards-compatibility loader in pecotmr. - twas_weights_conversion.R: convert legacy weight RDS (univariate and multicontext) into one pecotmr TwasWeights per gene. Weights and CV metrics are copied unchanged; carries traitPos (gene TSS) for cTWAS block placement; writes the xqtl_meta row and the context -> data type rows. - twas_ctwas.ipynb: new optional twas_weights_conversion workflow (skip it for weights already in TwasWeights format), with Input/Output/MWE docs and the Overview updated; MWE commands now use the converted example tables. - ctwas_1: --gwas_study now selects the GWAS (cTWAS takes one per LD block) and the run fails early if several remain; new --ctwas_method (weights with several methods and no ensemble); pass the chromosome's genotype prefix as --ld-sketch so ctwasPipeline sees a real LD file rather than a "<chrom-meta>" placeholder. - ctwas_assemble.R: drop study elements and LD blocks with no GWAS variants (they carry no genome build and broke combineGwasSumStats). Tested: conversion is exact for eQTL, pQTL and sQTL weights (weights compared element by element); TWAS runs on converted weights for 8 studies in the APOE region; cTWAS assemble passes. cTWAS parameter estimation and fine-mapping are not yet tested. Co-Authored-By: Claude Opus 5.5 <[email protected]> Claude-Session: https://claude.ai/code/session_01Csueazw8VDi9a3mhLvgjQg
- fine_mapping_conversion.R: convert legacy FunGen-xQTL QTL fine-mapping (*.univariate_bvsr.rds) into pecotmr QtlFineMappingResult, one per study and gene, with the same study/context labels as twas_weights_conversion. - twas_ctwas.ipynb: new optional fine_mapping_conversion step; twas passes --fine-mapping-result when --fine_mapping_dir has a file for the gene, plus --mr_method and --mr_pip_cutoff; Inf is passed to R as Inf; docs note that --rsq_cutoff 0 --rsq_pval_cutoff inf keeps every method.
twasWeightsRow fits -> methodFits; TwasWeights and QtlFineMappingResult study -> studyName.
gwas_sumstats_construct.R now exits cleanly without writing output when a block has no GWAS variants, and the ctwas step assembles from a manifest of the non-empty blocks only.
The twas step does not write per-gene .twas_weights.rds files, so ctwas_1 now defaults to the TwasWeights files listed in the TWAS manifest. Docs updated to match.
ctwas_finemap.R: add --missing-group-prior (default fill). Groups present in region data but without an estimated prior get the mean prior and variance of the estimated gene groups, with a warning; use stop to keep the old error. twas_ctwas.ipynb: cTWAS example commands use ld_meta_file.ctwas.tsv (20 blocks), since prior estimation needs at least two LD blocks.
New optional step in twas_ctwas.ipynb and gwas_sumstats_to_tsv.R: export per-LD-block pecotmr GwasSumStats RDS files of one study into one bgzipped, tabix-indexed TSV per chromosome plus the gwas_meta TSV used by --gwas_meta_data. Empty blocks are skipped; a variant present in two adjacent blocks is kept once, from the first block. The TWAS example also notes how to split a genome-wide scan across a cluster.
twas.R now takes comma-separated --twas-weights, --fine-mapping-result and --output lists for genes of the same LD block, computes the block LD from the reference panel once and subsets it per gene, and records a gene with no scorable (QTL, GWAS) pair as an empty result with a noResult attribute instead of stopping the block. The twas step calls twas.R once per region. Results are identical to the per-gene calls on the 19 genes of chr2:109502880-112720775; the block runs about 9x faster.
…ocks gwas_sumstats_to_tsv.R: when a variant appears in two GwasSumStats blocks, keep the copy that was measured (has an allele frequency) over one imputed by RAISS in that block, and fall back to the block that starts first. On AD_Bellenguez_main_2026 this changes the kept z for 964 of 16,826 shared variants.
Optional text file of block names whose copy of a variant shared with another block is kept first, ahead of the measured-over-imputed rule. Used to keep the larger added block of each overlapping pair in the RSS GWAS r4 block set; on AD_Bellenguez_main_2026 this changes 3,436 variants, all inside the listed overlap regions.
This branch has not been deployed
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
TWAS/cTWAS fixes for genome-wide runs, GwasSumStats input step, faster TWAS step
These changes came out of running the TWAS/cTWAS notebook genome-wide on the 2026 AD GWAS with the released FunGen-xQTL weights (pecotmr 0.9.3).
cTWAS
--twas_weightsis not given, use the TwasWeights files listed in the manifest as the default cTWAS weights.ctwas_finemap.R: new--missing-group-prior(defaultfill). A group that appears in the region data but got no prior estimate (too few regions carry its genes) now gets the mean prior and variance of the estimated gene groups, with a warning.--missing-group-prior stopkeeps the old error.ld_meta_file.ctwas.tsv(20 blocks), since prior estimation needs at least two LD blocks.New step
gwas_sumstats_to_tsvgwas_metafile for--gwas_meta_data.--prefer_blocksfile (used for the larger "added" block of each overlapping pair in the r4 block set); then the measured copy over a RAISS-imputed one (imputed variants have no allele frequency); then the block that starts first.Faster TWAS step
twas.Rnow accepts comma-separated--twas-weights,--fine-mapping-resultand--outputfor the genes of one LD block.twasstep calls it once per region instead of once per gene.noResultattribute instead of stopping the rest of its block.Testing
twas_weights_conversion,twas, and the threectwassteps with the 20-block LD meta) completes, and the fixture gene returns a TWAS result. The chr19 MWE fine-mapping also completes with filled priors.gwas_sumstats_to_tsvon the full AD_Bellenguez main / nobiobank / noproxy 2026 r4 output (1,358 blocks each). The main-study output has the same variants and z-scores as the existing per-chromosome 2026 GWAS, except one boundary variant whose sign differs between its two blocks.noResultoutput and don't stop their block.