diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 471239753..b86c3b725 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -36,6 +36,13 @@ jobs: export PIXI_HOME=${HOME}/.pixi curl -fsSL https://raw.githubusercontent.com/StatFunGen/pixi-setup/main/pixi-setup.sh | bash + # Companion retention change: StatFunGen/pecotmr#600. + # Use its immutable revision until a released pecotmr includes the fix. + - name: Install pecotmr with retained fSuSiE curves + run: >- + Rscript -e 'remotes::install_github("hsun3163/pecotmr@40b77e48a53d37b8b46fa897909585e6054a3c46", + upgrade = "never")' + - name: Script tests run: pytest tests/scripts -n logical --dist loadscope diff --git a/code/SoS/mnm_analysis/mnm_methods/mnm_regression.ipynb b/code/SoS/mnm_analysis/mnm_methods/mnm_regression.ipynb index e9dc2f16e..1c51deb95 100644 --- a/code/SoS/mnm_analysis/mnm_methods/mnm_regression.ipynb +++ b/code/SoS/mnm_analysis/mnm_methods/mnm_regression.ipynb @@ -1278,7 +1278,6 @@ { "cell_type": "code", "execution_count": null, - "id": "62f20d88", "metadata": { "kernel": "SoS" }, @@ -1302,6 +1301,8 @@ " --output ${_output}\n", "\n", "[fsusie_2]\n", + "# Empty preserves standalone behavior; Snakemake passes its chromosome scope.\n", + "parameter: chromosomes = []\n", "# Functional SuSiE per locus (fsusieR::susiF) over the pre-built QtlDataset.\n", "# fsusie-specific knobs ride on --method-args; --susie-top-pc>0 additionally\n", "# fine-maps each context's top principal components with univariate SuSiE\n", @@ -1321,13 +1322,16 @@ "import csv, json\n", "manifest = path(f\"{cwd:a}/fsusie/{name}.fsusie_region_manifest.tsv\")\n", "jobs = list(csv.DictReader(open(manifest), delimiter='\\t'))\n", + "if chromosomes:\n", + " allowed_chromosomes = {str(chrom).removeprefix(\"chr\") for chrom in chromosomes}\n", + " jobs = [job for job in jobs if str(job[\"chr\"]).removeprefix(\"chr\") in allowed_chromosomes]\n", "stop_if(len(jobs) == 0, \"fsusie: empty region manifest; check --region-name / --customized-association-windows.\")\n", "qtl_dataset = path(f\"{cwd:a}/qtl_dataset/{name}.qtl_dataset.rds\")\n", "fsusie_args = fsusie_method_args if fsusie_method_args else json.dumps({\"fsusie\": {\"prior\": prior, \"max_scale\": max_scale, \"min_purity\": min_purity, \"post_processing\": post_processing, \"max_SNP_EM\": max_SNP_EM, \"cor_small\": bool(small_sample_correction)}})\n", "input: qtl_dataset, for_each = \"jobs\"\n", - "output: f\"{cwd:a}/fsusie/{name}.{_jobs['region_id']}.fsusie.rds\"\n", - "task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f\"{step_name}_{_output:bn}\"\n", - "bash: expand = \"${ }\", stderr = f\"{_output:n}.stderr\", stdout = f\"{_output:n}.stdout\", container = container, entrypoint = entrypoint\n", + "output: f\"{cwd:a}/fsusie/{name}.{_jobs['region_id']}.fsusie.rds\", f\"{cwd:a}/fsusie/{name}.{_jobs['region_id']}.fsusie.exported.bed.gz\"\n", + "task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f\"{step_name}_{_output[0]:bn}\"\n", + "bash: expand = \"${ }\", stderr = f\"{_output[0]:n}.stderr\", stdout = f\"{_output[0]:n}.stdout\", container = container, entrypoint = entrypoint\n", " Rscript ${modular_script_dir}/pecotmr_integration/fine_mapping.R \\\n", " --qtl-dataset ${_input} \\\n", " --region ${_jobs['ld_block']} \\\n", @@ -1340,8 +1344,33 @@ " ${(\"--contexts \" + contexts) if contexts else \"\"} \\\n", " --method-args '${fsusie_args}' \\\n", " ${(\"--use-pca --n-pcs \" + str(susie_top_pc)) if susie_top_pc > 0 else \"\"} \\\n", - " --output ${_output}\n" - ] + " --output ${_output[0]}\n" + ], + "id": "62f20d88" + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "kernel": "SoS" + }, + "outputs": [], + "source": [ + "[fsusie_3]\n", + "# Assemble the same regional tables into one indexed context table.\n", + "input: output_from(\"fsusie_2\")\n", + "bed_inputs = paths([p for p in _input if str(p).endswith(\".exported.bed.gz\")])\n", + "output: f\"{cwd:a}/fsusie/{name}.exported.bed.gz\", f\"{cwd:a}/fsusie/{name}.exported.bed.gz.tbi\"\n", + "task: trunk_workers=1, trunk_size=1, mem=\"40G\", walltime=\"1h\", cores=1\n", + "bash: expand=\"${ }\", stderr=f\"{_output[0]}.stderr\", stdout=f\"{_output[0]}.stdout\", container=container, entrypoint=entrypoint\n", + " set -euo pipefail\n", + " {\n", + " gzip -cd ${bed_inputs[0]:q} | sed -n '1p'\n", + " gzip -cd ${bed_inputs:q} | awk '!/^#chr\\t/' | LC_ALL=C sort -k1,1V -k2,2n -k3,3n\n", + " } | bgzip -c > ${_output[0]:q}\n", + " tabix -f -p bed ${_output[0]:q}\n" + ], + "id": "fsusie-export" }, { "cell_type": "markdown", diff --git a/code/script/data_preprocessing/phenotype/prepare_qtl_manifest.py b/code/script/data_preprocessing/phenotype/prepare_qtl_manifest.py new file mode 100644 index 000000000..7e3ddc152 --- /dev/null +++ b/code/script/data_preprocessing/phenotype/prepare_qtl_manifest.py @@ -0,0 +1,77 @@ +#!/usr/bin/env python3 +"""Build the single-context phenotype manifest for pecotmr's QtlDataset. + +Every feature on the configured chromosomes points at the same phenotype BED: +``loadQtlDatasetFromManifest`` resolves one context to one phenotype matrix. +""" + +from __future__ import annotations + +import argparse +import csv +import gzip +import re +from pathlib import Path + + +def _open_text(path: Path): + if path.name.endswith((".gz", ".bgz")): + return gzip.open(path, "rt", newline="") + return path.open(newline="") + + +def _canonical(chromosome: str) -> str: + return re.sub(r"^chr", "", chromosome, flags=re.IGNORECASE).casefold() + + +def main() -> None: + parser = argparse.ArgumentParser( + description="Prepare a deterministic single-context QtlDataset manifest" + ) + parser.add_argument("--bed", required=True, type=Path) + parser.add_argument("--covariates", required=True, type=Path) + parser.add_argument("--context", required=True) + parser.add_argument("--chromosomes", required=True, nargs="+") + parser.add_argument("--phenotype-id-column", default="ID") + parser.add_argument("--phenotype-manifest", required=True, type=Path) + parser.add_argument("--region-ids", required=True, type=Path) + parser.add_argument("--sample-ids", required=True, type=Path) + args = parser.parse_args() + + bed = args.bed.resolve() + covariates = args.covariates.resolve() + order = {_canonical(c): i for i, c in enumerate(args.chromosomes)} + + rows = [] + with _open_text(bed) as handle: + reader = csv.reader(handle, delimiter="\t") + header = next(reader) + expected = ["#chr", "start", "end", args.phenotype_id_column] + if header[:4] != expected: + raise ValueError(f"{bed} must begin with " + ", ".join(expected)) + for chromosome, start, end, feature_id in (row[:4] for row in reader): + if _canonical(chromosome) in order: + rows.append([chromosome, int(start), int(end), feature_id, + str(bed), args.context, str(covariates)]) + if not rows: + raise ValueError(f"No features on the configured chromosomes: {bed}") + rows.sort(key=lambda row: (order[_canonical(row[0])], row[1], row[2], row[3])) + + with _open_text(covariates) as handle: + covariate_header = next(csv.reader(handle, delimiter="\t")) + samples = sorted(set(header[4:]) & set(covariate_header[1:])) + if not samples: + raise ValueError("Phenotype BED and covariate file share no sample IDs") + + for path in [args.phenotype_manifest, args.region_ids, args.sample_ids]: + path.parent.mkdir(parents=True, exist_ok=True) + with args.phenotype_manifest.open("w", newline="") as handle: + writer = csv.writer(handle, delimiter="\t", lineterminator="\n") + writer.writerow(["#chr", "start", "end", "ID", "path", "cond", "cov_path"]) + writer.writerows(rows) + args.region_ids.write_text("".join(f"{row[3]}\n" for row in rows)) + args.sample_ids.write_text("".join(f"{sample}\n" for sample in samples)) + + +if __name__ == "__main__": + main() diff --git a/code/script/pecotmr_integration/fine_mapping.R b/code/script/pecotmr_integration/fine_mapping.R index d1ef895e3..482b86414 100644 --- a/code/script/pecotmr_integration/fine_mapping.R +++ b/code/script/pecotmr_integration/fine_mapping.R @@ -60,6 +60,67 @@ suppressPackageStartupMessages({ library(jsonlite) }) +# The 19-column FunGen_xQTL_epi.bulk.exported.bed schema plus one packed +# grid_band_halfwidth column. Native fsusieR bands are curve +/- 3*sqrt(tt). +export_fsusie_bed <- function(x, path, region) { + columns <- c("#chr","start","end","variant_ID","event_ID","region_ID", + "maf","PIP","cs_coverage_0.95","TADB_start","TADB_end", + "grid_resolution","cs_id","cs_root","grid_positions", + "grid_effects","epi_mark_positions","epi_mark_names","epi_mark_effects", + "grid_band_halfwidth") + pieces <- list() + pack <- function(z) paste(z,collapse=";") + rg <- GenomicRanges::GRanges(region) + region_id <- paste0(as.character(GenomicRanges::seqnames(rg)),"_", + GenomicRanges::start(rg),"_",GenomicRanges::end(rg)) + for (i in which(x$method=="fsusie")) { + e <- x[i, ]; fit <- getSusieFit(e) + if (is.null(fit$fitted_func) || is.null(fit$cred_band) || is.null(fit$outing_grid)) + stop("Functional fields absent; rerun this region with curve retention") + tl <- as.data.frame(getTopLoci(e, signalCutoff = 0)) + grid <- as.numeric(fit$outing_grid) + positions <- as.numeric(fit$trait_positions) + probes <- as.character(fit$trait_names) + context <- as.character(x$context[i]) + stopifnot(length(grid)>1L,all(diff(grid)>0),length(positions)==length(probes)) + # Native fSuSiE CSs and curves use the same effect index. + for (l in seq_along(fit$cs)) { + members <- fit$cs[[l]] + if (!length(members)) next + vids <- if (is.numeric(members)) getVariantIds(e)[as.integer(members)] else as.character(members) + t <- tl[match(vids,tl$variant_id),,drop=FALSE] + stopifnot(!anyNA(t$variant_id)) + curve <- as.numeric(fit$fitted_func[[l]]) + band <- fit$cred_band[[l]] + stopifnot(length(curve)==length(grid),is.matrix(band), + identical(dim(band),c(2L,length(grid)))) + midpoint <- as.numeric((band[1,]+band[2,])/2) + stopifnot(isTRUE(all.equal(midpoint,curve,check.attributes=FALSE))) + halfwidth <- as.numeric((band[1,]-band[2,])/2) + # Linear interpolation and constant endpoint extrapolation match + # interpolate_effect_estimates in Update_new_epi_export.ipynb. + interpolated <- approx(grid,curve,xout=positions,rule=2,ties="ordered")$y + cs_id <- paste(context,region_id,l,sep=":") + pieces[[length(pieces)+1L]] <- data.frame( + `#chr`=sub("^chr","",t$chrom), start=t$pos-1L,end=t$pos, + variant_ID=t$variant_id,event_ID=context,region_ID=region_id, + maf=pmin(t$af,1-t$af),PIP=t$pip,cs_coverage_0.95=l, + TADB_start=GenomicRanges::start(rg),TADB_end=GenomicRanges::end(rg), + grid_resolution=length(grid),cs_id=cs_id,cs_root=cs_id, + grid_positions=pack(grid),grid_effects=pack(curve), + epi_mark_positions=pack(positions),epi_mark_names=pack(probes), + epi_mark_effects=pack(interpolated), + grid_band_halfwidth=pack(halfwidth),check.names=FALSE) + } + } + out <- if (length(pieces)) do.call(rbind,pieces) else + as.data.frame(setNames(replicate(length(columns),character(),simplify=FALSE),columns),check.names=FALSE) + out <- out[,columns,drop=FALSE] + if (nrow(out)) out <- out[order(out$`#chr`,out$start,out$end),,drop=FALSE] + f <- gzfile(path,"wt"); on.exit(close(f)) + write.table(out,f,sep="\t",quote=FALSE,row.names=FALSE,na="NA") +} + parser <- arg_parser("SuSiE fine-mapping over a pecotmr S4 input (QtlDataset or GwasSumStats)") parser <- add_argument(parser, "--qtl-dataset", help = "Path to a QtlDataset RDS (QTL mode)", @@ -418,7 +479,7 @@ if (has_gwas || has_qss) { label <- if (has_region) paste0("region '", argv$region, "'") else paste0("gene '", argv$gene_id, "'") run_fm <- function() if (has_region) { - do.call(fineMappingPipeline, c(qtl_args, list(region = argv$region))) + do.call(fineMappingPipeline, c(qtl_args, list(region = GenomicRanges::GRanges(argv$region)))) } else { do.call(fineMappingPipeline, c(qtl_args, list(traitId = argv$gene_id, @@ -475,3 +536,7 @@ if (has_gwas && !is.null(res)) { }, error = function(e) message("fine_mapping.R: GWAS summary skipped (", conditionMessage(e), ")"))) } + +if (has_qtl && has_region && "fsusie" %in% methods) { + export_fsusie_bed(res, sub("\\.rds$", ".exported.bed.gz", argv$output), argv$region) +} diff --git a/code/snakemake/README.md b/code/snakemake/README.md index 052d2c373..99e4327ad 100644 --- a/code/snakemake/README.md +++ b/code/snakemake/README.md @@ -280,3 +280,31 @@ Do not accidentally include those unless the user explicitly wants them. non-trivial test package is 288K. 4. If runtime verification is requested, run `run_mwe_xqtl_core.sh` against the external MWE data and run `run_nontrivial_tensorqtl_susie.sh`. + +### Functional SuSiE + +The shared `fsusie_finemapping` target prepares a single-context QtlDataset, +selects TAD regions, and calls the canonical `mnm_regression.ipynb` `fsusie` +workflow. Configure the existing `finemapping` section, for example: + +```yaml +finemapping: + fsusie: + tad_list: /path/to/TADs.bed + phenotype_per_tad: 16 + cis_window: 0 + susie_top_pc: 1 + post_processing: TI +``` + +The fSuSiE workflow writes each regional RDS and a matching exported BED, then +combines the regional tables into `fsusie/.exported.bed.gz` with a tabix +index. This requires a pecotmr version retaining fSuSiE functional summaries +and trait names/positions in trimmed fits. Older results lacking those fields +must be refitted. The table contains credible-set variants, grid effects, +interpolated probe effects, and `grid_band_halfwidth`. The latter is a packed +vector aligned with the grid: reconstruct the saved native bands as +`grid_effects -/+ grid_band_halfwidth`. These are native fSuSiE bands, not +standard errors. A region without a credible set contributes only its header. +TI is the default. `susie_top_pc` controls additional univariate PC fits; this +workflow does not invoke the TWAS-weight training steps. diff --git a/code/snakemake/Snakefile b/code/snakemake/Snakefile index f34b85503..7c3f1c102 100644 --- a/code/snakemake/Snakefile +++ b/code/snakemake/Snakefile @@ -436,6 +436,20 @@ def get_input_plink(wc): return f"{CWD}/data_preprocessing/genotype/{CONVERTED_PLINK_BASENAME}.bed" +def sos_mem_arg(mem_mb): + """Forward at least 32768 MB to a SoS task in the config's native unit.""" + return f"{max(32768, int(mem_mb))}M" + +def get_fsusie_region_list(theme): + fsusie_cfg = config["finemapping"].get("fsusie", {}) + return ( + f"{CWD}/finemapping/{theme}/fsusie/" + f"{Path(get_phenotype_bed_by_theme(theme)).name}." + f"{Path(fsusie_cfg.get('tad_list', '')).name}." + f"{fsusie_cfg.get('phenotype_per_tad', 2)}_pheno_per_region.region_list" + ) + + # ── Include script-backed rule modules ───────────────────────────────────────────── include: "rules/00_phenotype_preprocessing.smk" include: "rules/01_molecular_phenotypes.smk" @@ -444,6 +458,7 @@ include: "rules/03_sample_qc_pca.smk" include: "rules/04_phenotype_covariate_prep.smk" include: "rules/05_association_testing.smk" include: "rules/06_univariate_finemapping.smk" +include: "rules/09_fsusie_finemapping.smk" # ============================================================ @@ -554,3 +569,12 @@ rule finemapping: "{cwd}/finemapping/{theme}/susie_twas/.done_susie_twas", cwd=CWD, theme=THEMES, ), + + +rule fsusie_finemapping: + """Run functional SuSiE over TAD regions for every configured theme.""" + input: + expand( + "{cwd}/finemapping/{theme}/fsusie/.done_fsusie", + cwd=CWD, theme=THEMES, + ), diff --git a/code/snakemake/rules/09_fsusie_finemapping.smk b/code/snakemake/rules/09_fsusie_finemapping.smk new file mode 100644 index 000000000..d6fb84526 --- /dev/null +++ b/code/snakemake/rules/09_fsusie_finemapping.smk @@ -0,0 +1,180 @@ +# ============================================================ +# Rule Module 09: Functional Fine-mapping (fSuSiE) (script-backed) +# ============================================================ +# Covers: fSuSiE fine-mapping of functional phenotypes over TAD regions +# +# SoS notebooks called (script-backed wrappers in pipeline/): +# - phenotype_formatting.ipynb (phenotype_annotate_by_tad) +# - mnm_regression.ipynb (qtl_dataset_construct, fsusie) +# ============================================================ + +# ------------------------------------ +# Step 9.1 — Build the single-context QtlDataset inputs +# ------------------------------------ +rule fsusie_qtl_manifest: + """Describe one phenotype matrix, its features, and aligned samples for pecotmr.""" + input: + phenotype = lambda wc: get_phenotype_bed_by_theme(wc.theme), + hidden_factors = lambda wc: get_hidden_factors(wc), + output: + manifest = CWD + "/finemapping/{theme}/fsusie/qtl_dataset/{theme}.phenotype_manifest.tsv", + region_ids = CWD + "/finemapping/{theme}/fsusie/qtl_dataset/{theme}.region_ids.txt", + sample_ids = CWD + "/finemapping/{theme}/fsusie/qtl_dataset/{theme}.sample_ids.txt", + params: + script = MODULAR_SCRIPT_DIR + "/data_preprocessing/phenotype/prepare_qtl_manifest.py", + chromosomes = " ".join(config["chromosomes"]), + phenotype_id_column = lambda wc: _theme_cfg(wc.theme).get("phenotype_id_column", "ID"), + threads: 1 + resources: + mem_mb = config["resources"]["default"]["mem_mb"], + runtime = config["resources"]["default"]["runtime"], + shell: + """ + python3 {params.script} \ + --bed {input.phenotype} \ + --covariates {input.hidden_factors} \ + --context {wildcards.theme} \ + --chromosomes {params.chromosomes} \ + --phenotype-id-column {params.phenotype_id_column} \ + --phenotype-manifest {output.manifest} \ + --region-ids {output.region_ids} \ + --sample-ids {output.sample_ids} + """ + + +# ------------------------------------ +# Step 9.2 — Construct the QtlDataset once per context +# ------------------------------------ +rule fsusie_qtl_dataset: + """Build the QtlDataset shared by every regional fSuSiE task.""" + input: + genotype = lambda wc: config["finemapping"].get("genotype_file", get_plink_qc_bed()), + manifest = CWD + "/finemapping/{theme}/fsusie/qtl_dataset/{theme}.phenotype_manifest.tsv", + sample_ids = CWD + "/finemapping/{theme}/fsusie/qtl_dataset/{theme}.sample_ids.txt", + hidden_factors = lambda wc: get_hidden_factors(wc), + output: + qtl_dataset = CWD + "/finemapping/{theme}/fsusie/qtl_dataset/{theme}.qtl_dataset.rds", + params: + sos_bin = SOS_BIN, + sos_sched = sos_sched("qtl_dataset_construct"), + notebooks_dir = NOTEBOOKS, + modular_script_dir = MODULAR_SCRIPT_DIR, + outdir = CWD + "/finemapping/{theme}/fsusie", + maf = config["finemapping"]["maf"], + mac = config["association"]["mac_threshold"], + sos_mem = sos_mem_arg(config["resources"]["finemapping"]["mem_mb"]), + sos_walltime = sos_walltime_arg(config["resources"]["finemapping"]["runtime"], sos_queue("qtl_dataset_construct")), + dry_run = DRY_RUN_SOS, + threads: 1 + resources: + mem_mb = config["resources"]["finemapping"]["mem_mb"], + shell: + """ + {params.sos_bin} run {params.notebooks_dir}/mnm_regression.ipynb qtl_dataset_construct \ + --cwd {params.outdir} \ + --name {wildcards.theme} \ + --study {wildcards.theme} \ + --genoFile {input.genotype} \ + --phenoFile {input.manifest} \ + --covFile {input.hidden_factors} \ + --transpose-covariates \ + --maf-cutoff {params.maf} \ + --mac-cutoff {params.mac} \ + --keep-samples {input.sample_ids} \ + --mem {params.sos_mem} \ + --walltime {params.sos_walltime} \ + --modular-script-dir {params.modular_script_dir} \ + --numThreads {threads} {params.dry_run} {params.sos_sched} + """ + + +# ------------------------------------ +# Step 9.3 — TAD regions holding enough phenotype features +# ------------------------------------ +rule fsusie_region_list: + """Keep the TADs that contain at least phenotype_per_tad phenotype features.""" + input: + phenotype = lambda wc: get_phenotype_bed_by_theme(wc.theme), + tad_list = config["finemapping"].get("fsusie", {}).get("tad_list", "") or [], + output: + region_list = CWD + "/finemapping/{theme}/fsusie/{region_list_base}.region_list", + params: + sos_bin = SOS_BIN, + notebooks_dir = NOTEBOOKS, + modular_script_dir = MODULAR_SCRIPT_DIR, + outdir = CWD + "/finemapping/{theme}/fsusie", + phenotype_per_tad = config["finemapping"].get("fsusie", {}).get("phenotype_per_tad", 2), + dry_run = DRY_RUN_SOS, + threads: 1 + resources: + mem_mb = config["resources"]["default"]["mem_mb"], + shell: + """ + {params.sos_bin} run {params.notebooks_dir}/phenotype_formatting.ipynb phenotype_annotate_by_tad \ + --cwd {params.outdir} \ + --phenoFile {input.phenotype} \ + --TAD-list {input.tad_list} \ + --phenotype-per-tad {params.phenotype_per_tad} \ + --modular-script-dir {params.modular_script_dir} \ + --numThreads {threads} {params.dry_run} + """ + + +# ------------------------------------ +# Step 9.4 — fSuSiE per TAD region over the QtlDataset +# ------------------------------------ +rule fsusie: + """Fine-map functional phenotypes across each retained TAD region.""" + input: + genotype = lambda wc: config["finemapping"].get("genotype_file", get_plink_qc_bed()), + manifest = CWD + "/finemapping/{theme}/fsusie/qtl_dataset/{theme}.phenotype_manifest.tsv", + qtl_dataset = CWD + "/finemapping/{theme}/fsusie/qtl_dataset/{theme}.qtl_dataset.rds", + region_list = lambda wc: get_fsusie_region_list(wc.theme), + hidden_factors = lambda wc: get_hidden_factors(wc), + output: + done = CWD + "/finemapping/{theme}/fsusie/.done_fsusie", + table = CWD + "/finemapping/{theme}/fsusie/fsusie/{theme}.exported.bed.gz", + index = CWD + "/finemapping/{theme}/fsusie/fsusie/{theme}.exported.bed.gz.tbi", + params: + sos_bin = SOS_BIN, + sos_sched = sos_sched("fsusie"), + notebooks_dir = NOTEBOOKS, + modular_script_dir = MODULAR_SCRIPT_DIR, + outdir = CWD + "/finemapping/{theme}/fsusie", + pip_cutoff = config["finemapping"]["pip_cutoff"], + cis_window = config["finemapping"].get("fsusie", {}).get("cis_window", 0), + post_processing = config["finemapping"].get("fsusie", {}).get("post_processing", "TI"), + susie_top_pc = config["finemapping"].get("fsusie", {}).get("susie_top_pc", 0), + chromosomes = " ".join(config["chromosomes"]), + coverage = " ".join(str(x) for x in config["finemapping"]["coverage"]), + seed = config.get("analysis", {}).get("seed", 999), + sos_mem = sos_mem_arg(config["resources"]["finemapping"]["mem_mb"]), + sos_walltime = sos_walltime_arg(config["resources"]["finemapping"]["runtime"], sos_queue("fsusie")), + dry_run = DRY_RUN_SOS, + threads: 1 + resources: + mem_mb = config["resources"]["finemapping"]["mem_mb"], + shell: + """ + {params.sos_bin} run {params.notebooks_dir}/mnm_regression.ipynb fsusie \ + --cwd {params.outdir} \ + --name {wildcards.theme} \ + --genoFile {input.genotype} \ + --phenoFile {input.manifest} \ + --covFile {input.hidden_factors} \ + --customized-association-windows {input.region_list} \ + --cis-window {params.cis_window} \ + --chromosomes {params.chromosomes} \ + --susie-top-pc {params.susie_top_pc} \ + --post-processing {params.post_processing} \ + --pip-cutoff {params.pip_cutoff} \ + --coverage {params.coverage} \ + --seed {params.seed} \ + --mem {params.sos_mem} \ + --walltime {params.sos_walltime} \ + --modular-script-dir {params.modular_script_dir} \ + --numThreads {threads} {params.dry_run} {params.sos_sched} + status=$? + if [ "$status" -ne 0 ]; then exit "$status"; fi + touch {output.done} + """ diff --git a/tests/notebooks/mnm_analysis/mnm_methods/test_mnm_regression.py b/tests/notebooks/mnm_analysis/mnm_methods/test_mnm_regression.py index bb6a5b4c9..ad8774aac 100644 --- a/tests/notebooks/mnm_analysis/mnm_methods/test_mnm_regression.py +++ b/tests/notebooks/mnm_analysis/mnm_methods/test_mnm_regression.py @@ -83,3 +83,53 @@ def test_mnm(run_sos, read_rds, repo_root, qtl_mini, tmp_path): exp = repo_root / "tests/fixtures/mnm_regression/expected" assert_matches_expected(fmr, exp / "multicontext_bvsr.rds", mode="tolerant", rtol=1e-6, atol=1e-8) + + +def test_fsusie_ti_export(run_sos, read_rds, repo_root, qtl_mini, tmp_path): + """The native fSuSiE workflow writes one indexed table with curves and bands.""" + import csv + import gzip + import json + import subprocess + + cwd = tmp_path / "fsusie" + windows = tmp_path / "functional_region.bed" + windows.write_text("#chr\tstart\tend\tID\nchr22\t10000000\t18000000\tfunctional_region\n") + p = run_sos( + repo_root / "pipeline/mnm_regression.ipynb", + "qtl_dataset_construct+fsusie", + { + "name": "test_study", "cwd": cwd, + "genoFile": qtl_mini / "protocol_example.genotype.chr22.bed", + "phenoFile": qtl_mini / "protocol_example.pheno_manifest_context.tsv", + "covFile": qtl_mini / "example_covariates.tsv", + "customized-association-windows": windows, + "transpose-covariates": True, "seed": 1, + "susie-top-pc": 1, "mem": "40G", + "fsusie-method-args": json.dumps({"fsusie": { + "post_processing": "TI", "max_scale": 4, "L": 2, + "max_SNP_EM": 20, "verbose": False}}), + "modular_script_dir": repo_root / "code/script", + }, cwd=repo_root, timeout=1800) + assert p.returncode == 0, p.stdout + p.stderr + fits = list((cwd / "fsusie").glob("*.fsusie.rds")) + assert len(fits) == 1 + info = read_rds(fits[0]) + assert info["class"] == "QtlFineMappingResult" + assert info["nrow"] == 4 # one joint fit + one PC per context + table = cwd / "fsusie/test_study.exported.bed.gz" + assert table.with_suffix(table.suffix + ".tbi").exists() + with gzip.open(table, "rt") as handle: + reader = csv.DictReader(handle, delimiter="\t") + assert len(reader.fieldnames) == 20 + rows = list(reader) + assert rows, "Fixture must exercise a nonempty credible-set export" + for row in rows: + n = int(row["grid_resolution"]) + assert len(row["grid_effects"].split(";")) == n + widths = [float(v) for v in row["grid_band_halfwidth"].split(";")] + assert len(widths) == n and min(widths) >= 0 and max(widths) > 0 + assert len(row["epi_mark_names"].split(";")) == len(row["epi_mark_effects"].split(";")) + p = subprocess.run(["tabix", str(table), "22:10000000-18000000"], + capture_output=True, text=True) + assert p.returncode == 0 and p.stdout.strip()