Fix MWE bugs across genotype preprocessing, enrichment, EMS, TWAS/colocalization, and molecular phenotype notebooks#1387
Merged
Conversation
added 12 commits
July 8, 2026 03:32
…uired cohort/context parameters GWAS_QC, PCA, genotype_formatting, genotype_alignment, and TensorQTL had a 'parameter: name = ""' (or similarly unused) placeholder that was never wired to the documented --name example, so the MWE commands silently used an empty run tag. Replaced with explicit required parameters (--cohort, and --cohort/--context/--modality for TensorQTL) that are actually used to build the output name, matching the protocol_example naming convention.
…essing function name eoo_enrichment's documented example passed --name, which the [global] step never declared (it only had parameter: name = 'eoo'); replaced with the real required parameters --trait/--annotation-name so the example runs as written. sldsc_enrichment's docs/text referenced pecotmr::sldsc_postprocessing_pipeline, which does not exist in the installed pecotmr version; the real exported function is sldscPostprocessingPipeline. Updated all references. Note: the postprocess/meta_subset steps still fail on a deeper argument-signature mismatch (old flat-argument calling convention vs. the new S4-object-based SldscData API) which is not fixed here and needs a follow-up.
… prediction model-path default - Default --config-dir pointed at code/xqtl_modifier_score, which does not exist; the real location is code/SoS/xqtl_modifier_score. - The [train] bash step ran 'touch $[_output]' unconditionally after the python call with no 'set -e', so a failed training run was still reported as successful. Added set -e to the block. - data_config.yaml used ../.. for paths relative to the script, but the script actually runs from code/SoS/xqtl_modifier_score (3 levels deep, not 2), so every input path resolved incorrectly. Fixed to ../../... - ems_prediction's documented --model-path pointed at output/ems_training/protocol_example_chr2_scEEMS_model.joblib, which is never produced; corrected to the real path under output/xqtl_modifier_score/protocol_example/model_results/. Validated end-to-end: train produces a real cross-validated model (joblib), predict produces real per-variant EMS scores.
…fix ld_prune_reference variable name rss_ld_sketch's documented --cohort-id protocol_example. had a stray trailing dot, producing double-dot output filenames (protocol_example..pgen etc). Removed the trailing dot from both occurrences and regenerated the chr22 outputs to confirm clean single-dot filenames. ld_prune_reference referenced an undefined variable 'genotype' instead of the actual parameter 'genotype_list' when reading the genotype file list.
…ion notebooks
- mnm_regression: drop unavailable susieInf method from default twas_methods.
- mnm_postprocessing: rename unused 'name' placeholder to required --study,
actually wired into the output name.
- univariate_fine_mapping_twas_vignette: documented command skipped the
required qtl_dataset_construct step and pointed at a phenotype file that
does not match the multi-trait manifest schema; fixed to chain
qtl_dataset_construct+susie_twas with the correct manifest-schema files.
- SuSiE_enloc: fixed analysis_name_prefix matching, which used
split('.')[0] assuming the wrong filename format; corrected to a
case-insensitive substring match against the real
{prefix}.{cohort}.{analysis_name}.{gene}...rds naming; also corrected the
documented enrichment output structure/example to match real output.
- twas_ctwas: added a fallback that writes explicit empty placeholder
outputs when a region produces none, so downstream steps do not break on
missing files.
- legacy_twas_weights_convert.R: fixed output field name cvPerformance ->
cvResult to match the consumer's expected schema.
…nfo in pseudobulk BED construction to_midpoint_bed assumed every peak ID encodes chr-start-end and could be split on '-', which fails for gene-ID-based peaks (e.g. plain Ensembl IDs). Added a GTF-based gene-position lookup (load_gene_positions) used as a fallback to build the BED record (TSS-based single-base coordinate) when the peak ID isn't itself a coordinate triplet. snRNAseq_preprocessing: no functional change, only re-saved with escaped unicode in comments (box-drawing separators / arrows).
…otmr API renames; parameter placeholders mash_preprocessing: load_multitrait_R_sumstat requires extract_top_loci, which is not present in the installed pecotmr version. Added a direct, conceptually-equivalent construction of mash_input from the toy sumstats_db.rds (per-condition z-scores merged on common variants) so the step can run against the MWE data. Also renamed merge_mash_data -> mergeMashData and filter_invalid_summary_stat -> filterInvalidSummaryStat to match the installed pecotmr API. MRAID_QTL, fastenloc_dap, polyfun (graveyard/deprecated notebooks): replaced placeholder 'parameter: name = "demo"' with the real required parameters actually used to build the run name, and fixed a 'paramter' typo.
Per feedback, simplified back from the --cohort/--context/--modality split to a single --name parameter across TensorQTL, GWAS_QC, PCA, genotype_formatting, and genotype_alignment. Kept it as a required parameter (no empty-string default) instead of the original 'parameter: name = ""', so it can no longer be silently left blank/unused -- it must be explicitly supplied, restoring the simpler single-flag CLI.
…n) untouched Per feedback, these deprecated/graveyard notebooks should not be modified. Reverted them back to their state prior to the placeholder parameter fix; mash_preprocessing.ipynb (not a graveyard notebook) keeps its fix.
…tead of mismatched 8-condition computed prior
…bsolutize posterior_vhat_files - posterior_1: posterior_vhat_files was not converted to an absolute path before the R block chdir's into --cwd, so the documented relative --posterior-vhat-files example silently failed with "cannot open compressed file". Now absolutized like mash_model already was. - mash_posterior_contrast doc example was invoking the wrong workflow target/parameters entirely (--analysis-units/--mash-model/--posterior-input, which belong to the posterior step) and referenced a --posterior-input flag that does not exist on this target. Replaced with the correct two-step posterior -> mash_posterior_contrast invocation using --posterior-file and --sum-file (tab-separated id/path manifests), matching what the step actually implements. Verified end-to-end with a clean run producing real contrast output (posterior_sum.csv/png).
…g compatibility bugs - Add parameter: dapars_path to [UTR_reference] and [APAmain] steps in apa_calling.ipynb, defaulting to the real in-repo location code/SoS/molecular_phenotypes/calling/apa. This replaces a hardcoded absolute path to a different author's machine (/mnt/mfs/statgen/ls3751/github/xqtl-protocol/code/Dapars2_Multi_Sample.py) which was broken for anyone else running the pipeline. - Invoke gtf2bed12.py/DaPars_Extract_Anno.py/Dapars2_Multi_Sample.py via 'python2 <dapars_path>/<script>' instead of bare command names, since these scripts are not on PATH and are not marked executable. - Fix documented example commands (apa.ipynb, apa_calling.ipynb) which referenced code/molecular_phenotypes/... (missing the SoS/ path segment) for --dapars-path. - Fix gtf2bed12.py: the gene-type regex only matched GENCODE-style gene_type="..." attributes, silently producing 0 output on the toy chr22.gtf which uses Ensembl-style gene_biotype="...". Broadened the regex to gene_(?:bio)?type to accept both, backward compatible. Verified: UTR_reference now extracts 3766 3'UTR regions (was 0). - Add --chr-prefix true to the documented APAmain example: the toy wig files already contain 'chr'-prefixed chromosome names, but the default chr_prefix=False tells Dapars2_Multi_Sample.py to add a 'chr' prefix itself, producing 'chrchr22' which never matches and silently drops all coverage (visible as 'no wig: ...' warnings + KeyError in worker processes). Verified via clean end-to-end re-run of UTR_reference and APAmain with the corrected documented commands: no path/permission/parsing errors, real 3'UTR regions extracted, and APAmain completes without warnings or KeyErrors and writes a well-formed result file. Confirmed (via a temporary, reverted threshold override) that the toy dataset's read depth is simply below the pipeline's default Coverage_threshold=10 / least_pass_coverage_percentage=0.3, which is why the final PDUI table has 0 passing genes with default thresholds on this small toy subsample -- this is a data-scale/threshold limitation of the MWE, not a code bug. Also installed scipy for the system python2 via 'pip2 install --user scipy' (user-scoped), required by Dapars2_Multi_Sample.py.
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
2 participants
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.
This PR fixes a batch of genuine bugs found while systematically re-running the documented minimal working examples (MWEs) for each pipeline notebook. All fixes are targeted, code-level corrections validated by actually re-running the affected commands against the toy/protocol example data.