StatFunGen · gaow · Jul 21, 2026 · Jul 13, 2026 · Jul 14, 2026 · Jul 14, 2026
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
@@ -0,0 +1,34 @@
+name: Continuous Integration
+
+on:
+  pull_request:
+    paths:
+      - 'code/script/**'
+      - 'code/SoS/**'
+      - 'tests/**'
+  workflow_dispatch:
+
+concurrency:
+  group: wrapper-tests-${{ github.ref }}
+  cancel-in-progress: true
+
+jobs:
+  ci:
+    name: CI (${{ matrix.runner }})
+    runs-on: ${{ matrix.runner }}
+    strategy:
+      fail-fast: false
+      matrix:
+        runner: [ubuntu-latest, macos-latest]
+    steps:
+      - name: Checkout repository
+        uses: actions/checkout@v7
+
+      - name: Set up pixi environment
+        uses: prefix-dev/setup-pixi@v0.10.0
+
+      - name: Script tests
+        run: pixi run test_scripts
+
+      - name: Notebook tests
+        run: pixi run test_notebooks
diff --git a/.gitignore b/.gitignore
@@ -55,3 +55,5 @@ renovated_code/snakemake/archive/
 renovated_code/snakemake/dryrun/bin/
 
 dev/
+
+pixi.lock
diff --git a/code/SoS/association_scan/qtl_association_postprocessing.ipynb b/code/SoS/association_scan/qtl_association_postprocessing.ipynb
@@ -155,44 +155,30 @@
    "source": [
     "[global]\n",
     "parameter: cwd = path(\".\")\n",
-    "#path to pecotmr repo\n",
-    "parameter: pecotmr_path = path\n",
-    "parameter: gene_coordinates = path\n",
+    "# code/script dir that holds the pecotmr_integration wrappers\n",
+    "parameter: modular_script_dir = path\n",
     "parameter: output_dir = path\n",
-    "parameter: archive_dir = path\n",
-    "parameter: sub_dir = path\n",
+    "parameter: sub_dir = path(\".\")\n",
+    "# labels for the assembled QtlSumStats\n",
+    "parameter: study = \"study\"\n",
+    "parameter: context = \"context\"\n",
+    "parameter: genome = \"hg38\"\n",
     "parameter: maf_cutoff = 0.01\n",
     "parameter: cis_window = 1000000\n",
-    "parameter: tss_dist_col = \"start_distance\"\n",
-    "parameter: tes_dist_col = \"end_distance\"\n",
     "parameter: af_col = \"af\"\n",
+    "parameter: pvalue_col = \"pvalue\"\n",
     "parameter: molecular_id_col = \"molecular_trait_object_id\"\n",
-    "# This is for selecting the subset of data to process on protential signals \n",
-    "# assuming we drop those above this threshold\n",
-    "# This might lead to underestimates in qvalue method since qvalue < 0.05 may contain pvalue > 0.05\n",
+    "# drop pairs with p above this before assembling (keeps the object small; the\n",
+    "# per-gene min variant is always retained so the Bonferroni min is exact)\n",
     "parameter: pvalue_cutoff = 0.05\n",
-    "# This is used for both event and variant level significance filter\n",
     "parameter: fdr_threshold = 0.05\n",
-    "# eg \"pvalue\" for pvalue, \"pvalue_.*_interaction\" for interaction\n",
-    "parameter: pvalue_pattern = \"^pvalue$\"\n",
-    "# eg \"qvalue\" for pvalue, \"qvalue_interaction\" for interaction\n",
-    "parameter: qvalue_pattern = \"^qvalue$\"\n",
-    "# eg \"*.cis_qtl.regional.tsv.gz$\" for genetic effect via our pipeline TensorQTL.ipynb\n",
-    "# \"*.cis_qtl_top_assoc.txt.gz$\" for interaction genetic effect\n",
-    "parameter: regional_pattern = \"NULL\"\n",
-    "parameter: qtl_pattern = \"*.cis_qtl.pairs.tsv.gz$\"\n",
-    "parameter: n_variants_suffix = \"cis_n_variants_stats.tsv.gz\"\n",
+    "# list.files() regex patterns to locate the inputs in the work dir\n",
+    "parameter: regional_pattern = \"cis_qtl[.]regional[.]tsv[.]gz$\"\n",
+    "parameter: qtl_pattern = \"cis_qtl[.]pairs[.]tsv[.]gz$\"\n",
+    "parameter: n_variants_suffix = \"cis_n_variants_stats[.]tsv[.]gz$\"\n",
     "parameter: enable_archive = False\n",
-    "parameter: additional_pvalue_cols = \"\"\n",
     "\n",
-    "work_dir = f\"{cwd:a}/{sub_dir}\"\n",
-    "if sub_dir == path(\".\"):\n",
-    "    output_dir = f\"{output_dir:a}/{cwd:b}\"\n",
-    "    archive_dir = f\"{archive_dir:a}/{cwd:b}\"\n",
-    "    work_dir = f\"{cwd:a}\"\n",
-    "else:\n",
-    "    output_dir = f\"{output_dir:a}/{cwd:b}/{sub_dir}\"\n",
-    "    archive_dir = f\"{archive_dir:a}/{cwd:b}/{sub_dir}\""
+    "work_dir = cwd if str(sub_dir) == \".\" else path(f\"{cwd:a}/{sub_dir}\")\n"
    ]
   },
   {
@@ -204,59 +190,23 @@
    "outputs": [],
    "source": [
     "[default]\n",
-    "output: f\"{output_dir}/{cwd:b}_multiple_testing_consolidated.rds\"\n",
-    "task: trunk_workers = 1, tags = f'tensorqtl_postprocessing_{_output:n}'\n",
-    "R: expand = \"${ }\"\n",
-    "\n",
-    "    params <- list()\n",
-    "    params$workdir           <- \"${work_dir}\"\n",
-    "    params$maf_cutoff        <- ${maf_cutoff}\n",
-    "    params$cis_window        <- ${cis_window}\n",
-    "    params$pvalue_cutoff     <- ${pvalue_cutoff}\n",
-    "    params$fdr_threshold     <- ${fdr_threshold}\n",
-    "    params$gene_coordinates  <- \"${gene_coordinates:a}\"\n",
-    "    params$output_dir        <- \"${output_dir}\"\n",
-    "    params$archive_dir       <- \"${archive_dir}\"\n",
-    "    params$regional_pattern  <- \"${regional_pattern}\"\n",
-    "    params$n_variants_suffix <- \"${n_variants_suffix}\"\n",
-    "    params$qtl_pattern       <- \"${qtl_pattern}\"\n",
-    "    params$pvalue_pattern    <- \"${pvalue_pattern}\"\n",
-    "    params$qvalue_pattern    <- \"${qvalue_pattern}\"\n",
-    "    params$start_distance_col <- \"${tss_dist_col}\"\n",
-    "    params$end_distance_col   <- \"${tes_dist_col}\"\n",
-    "    params$af_col <- \"${af_col}\"\n",
-    "    params$molecular_id_col <- \"${molecular_id_col}\"\n",
-    "    enable_archive_val <- \"${enable_archive}\"\n",
-    "    if (enable_archive_val %in% c(\"True\", \"TRUE\", \"true\")) {\n",
-    "      params$enable_archive <- TRUE\n",
-    "    } else {\n",
-    "      params$enable_archive <- FALSE\n",
-    "    }\n",
-    "    message(sprintf(\"Archive setting - input: '%s', converted: %s\", enable_archive_val, params$enable_archive))    \n",
-    "    params$additional_pvalue_cols <- \"${additional_pvalue_cols}\"    \n",
-    "    convert_null_strings <- function(params) {\n",
-    "      if (is.list(params)) {\n",
-    "        # Apply the function to each element in the list\n",
-    "        result <- lapply(params, convert_null_strings)\n",
-    "        return(result)\n",
-    "      } else {\n",
-    "        # For non-list elements, check if it's the string \"NULL\"\n",
-    "        if (is.character(params) && params == \"NULL\") {\n",
-    "          return(NULL)\n",
-    "        } else {\n",
-    "          return(params)\n",
-    "        }\n",
-    "      }\n",
-    "    }\n",
-    "\n",
-    "    params <- convert_null_strings(params)\n",
-    "    source(\"${pecotmr_path}/inst/code/tensorqtl_postprocessor.R\")\n",
-    "    results <- hierarchical_multiple_testing_correction(params)\n",
-    "    write_results(results, params$output_dir, params$workdir, to_cwd = \"regional\")\n",
-    "    if (params$enable_archive) {\n",
-    "      archive_files(params)\n",
-    "    }\n",
-    "    saveRDS(results, ${_output:r})"
+    "# Hierarchical multiple-testing correction of the cis-QTL association results,\n",
+    "# via pecotmr::qtlAssociationPostprocess (the wrapper reads the per-gene regional\n",
+    "# + per-variant pairs, assembles a per-gene QtlSumStats, and writes the enriched\n",
+    "# consolidated RDS + per-method regional / significant-QTL / event / summary\n",
+    "# tables). Replaces the old source(pecotmr/inst/code/tensorqtl_postprocessor.R).\n",
+    "output: f\"{output_dir:a}/{cwd:b}.qtl_association_postprocessing.rds\"\n",
+    "bash: expand = \"${ }\"\n",
+    "    Rscript ${modular_script_dir:a}/pecotmr_integration/qtl_association_postprocessing.R \\\n",
+    "        --cwd ${work_dir:a} \\\n",
+    "        --regional-pattern \"${regional_pattern}\" \\\n",
+    "        --qtl-pattern \"${qtl_pattern}\" \\\n",
+    "        --n-variants-suffix \"${n_variants_suffix}\" \\\n",
+    "        --maf-cutoff ${maf_cutoff} --cis-window ${cis_window} \\\n",
+    "        --fdr-threshold ${fdr_threshold} --pvalue-cutoff ${pvalue_cutoff} \\\n",
+    "        --pvalue-col \"${pvalue_col}\" --af-col \"${af_col}\" \\\n",
+    "        --study \"${study}\" --context \"${context}\" --genome \"${genome}\" \\\n",
+    "        --output ${_output:a} --output-dir ${output_dir:a}\n"
    ]
   }
  ],

diff --git a/code/SoS/mnm_analysis/mnm_methods/mnm_regression.ipynb b/code/SoS/mnm_analysis/mnm_methods/mnm_regression.ipynb
@@ -313,7 +313,7 @@
    "source": [
     "### Multivariate analysis: mvSuSiE and mr.mash (`mnm`)\n",
     "\n",
-    "Fine-maps multiple molecular phenotypes jointly across a shared window using mvSuSiE with mr.mash priors. This step needs **multiple contexts** (e.g. the same genes measured under two or more conditions). The toy set ships a single context, so for this MWE we add a **second, synthetic context** (`example/colocboost_ctx2.bed.gz`, derived from context 1 with added noise — *demo data only, not a real measurement*) and a multi-context phenotype manifest. With that in place the `mnm` step now runs end-to-end on the toy data and produces `multicontext_bvsr.rds`, `multicontext_twas_weights.rds`, and `multicontext_data.rds` per gene."
+    "Fine-maps multiple molecular phenotypes jointly across a shared window using mvSuSiE with mr.mash priors. This step needs **multiple contexts** (e.g. the same genes measured under two or more conditions). The toy set ships a single context, so for this MWE we add a **second, synthetic context** (`example/colocboost_ctx2.bed.gz`, derived from context 1 with added noise \u2014 *demo data only, not a real measurement*) and a multi-context phenotype manifest. With that in place the `mnm` step now runs end-to-end on the toy data and produces `multicontext_bvsr.rds`, `multicontext_twas_weights.rds`, and `multicontext_data.rds` per gene."
    ]
   },
   {
@@ -578,60 +578,7 @@
     "    raise ValueError(f\"``skip_analysis_pip_cutoff`` should have either length 1 or length the same as phenotype files ({len(phenoFile)} in this case)\")\n",
     "\n",
     "# make it into an R List string\n",
-    "skip_analysis_pip_cutoff = [f\"'{y}'={x}\" for x,y in zip(skip_analysis_pip_cutoff, phenotype_names)]\n",
-    "    \n",
-    "def group_by_region(lst, partition):\n",
-    "    # from itertools import accumulate\n",
-    "    # partition = [len(x) for x in partition]\n",
-    "    # Compute the cumulative sums once\n",
-    "    # cumsum_vector = list(accumulate(partition))\n",
-    "    # Use slicing based on the cumulative sums\n",
-    "    # return [lst[(cumsum_vector[i-1] if i > 0 else 0):cumsum_vector[i]] for i in range(len(partition))]\n",
-    "    return partition\n",
-    "\n",
-    "import os\n",
-    "import pandas as pd\n",
-    "\n",
-    "def adapt_file_path(file_path, reference_file):\n",
-    "    \"\"\"\n",
-    "    Adapt a single file path based on its existence and a reference file's path.\n",
-    "\n",
-    "    Args:\n",
-    "    - file_path (str): The file path to adapt.\n",
-    "    - reference_file (str): File path to use as a reference for adaptation.\n",
-    "\n",
-    "    Returns:\n",
-    "    - str: Adapted file path.\n",
-    "\n",
-    "    Raises:\n",
-    "    - FileNotFoundError: If no valid file path is found.\n",
-    "    \"\"\"\n",
-    "    reference_path = os.path.dirname(reference_file)\n",
-    "\n",
-    "    # Check if the file exists\n",
-    "    if os.path.isfile(file_path):\n",
-    "        return file_path\n",
-    "\n",
-    "    # Check file name without path\n",
-    "    file_name = os.path.basename(file_path)\n",
-    "    if os.path.isfile(file_name):\n",
-    "        return file_name\n",
-    "\n",
-    "    # Check file name in reference file's directory\n",
-    "    file_in_ref_dir = os.path.join(reference_path, file_name)\n",
-    "    if os.path.isfile(file_in_ref_dir):\n",
-    "        return file_in_ref_dir\n",
-    "\n",
-    "    # Check original file path prefixed with reference file's directory\n",
-    "    file_prefixed = os.path.join(reference_path, file_path)\n",
-    "    if os.path.isfile(file_prefixed):\n",
-    "        return file_prefixed\n",
-    "\n",
-    "    # If all checks fail, raise an error\n",
-    "    raise FileNotFoundError(f\"No valid path found for file: {file_path}\")\n",
-    "\n",
-    "def adapt_file_path_all(df, column_name, reference_file):\n",
-    "    return df[column_name].apply(lambda x: adapt_file_path(x, reference_file))"
+    "skip_analysis_pip_cutoff = [f\"'{y}'={x}\" for x,y in zip(skip_analysis_pip_cutoff, phenotype_names)]\n"
    ]
   },
   {
@@ -718,8 +665,10 @@
     "# Fine-mapping SER pre-screen (0 = off, <0 = adaptive 3/nVariants).\n",
     "parameter: pip_cutoff_to_skip = 0.0\n",
     "# TWAS weight-learning knobs.\n",
-    "parameter: min_twas_maf = 0.01\n",
-    "parameter: min_twas_xvar = 0.01\n",
+    "# TWAS learns weights applied out-of-sample; keep an explicit MAF / X-variance\n",
+    "# cutoff for the TWAS pass (passed via the shared --maf-cutoff/--xvar-cutoff).\n",
+    "parameter: twas_maf_cutoff = 0.01\n",
+    "parameter: twas_xvar_cutoff = 0.01\n",
     "parameter: max_cv_variants = 5000\n",
     "parameter: cv_folds = 5\n",
     "parameter: cv_threads = 1\n",
@@ -743,15 +692,25 @@
     "qtl_dataset = path(f\"{cwd:a}/qtl_dataset/{name}.qtl_dataset.rds\")\n",
     "# Assemble the shared opt-in worker flags once (command construction only).\n",
     "seed_arg = f\"--seed {seed}\" if seed >= 0 else \"\"\n",
+    "# Non-MAF/xvar filters are shared by both passes (empty = use the dataset slots).\n",
     "filter_overrides = \" \".join(filter(None, [\n",
-    "    f\"--maf-cutoff {maf_cutoff}\"       if maf_cutoff   >= 0 else \"\",\n",
     "    f\"--mac-cutoff {mac_cutoff}\"       if mac_cutoff   >= 0 else \"\",\n",
-    "    f\"--xvar-cutoff {xvar_cutoff}\"     if xvar_cutoff  >= 0 else \"\",\n",
     "    f\"--imiss-cutoff {imiss_cutoff}\"   if imiss_cutoff >= 0 else \"\",\n",
     "    \"--drop-indel\"                     if drop_indel        else \"\",\n",
     "    f\"--keep-samples {keep_samples}\"   if keep_samples  != \".\" else \"\",\n",
     "    f\"--keep-variants {keep_variants}\" if keep_variants != \".\" else \"\",\n",
     "]))\n",
+    "# MAF / X-variance ride the shared --maf-cutoff/--xvar-cutoff, handled per pass:\n",
+    "# fine-mapping opt-in only (default off = dataset slot); TWAS always explicit\n",
+    "# (its portability cutoff) unless a general override is set, which wins for both.\n",
+    "fm_maf_xvar = \" \".join(filter(None, [\n",
+    "    f\"--maf-cutoff {maf_cutoff}\"   if maf_cutoff  >= 0 else \"\",\n",
+    "    f\"--xvar-cutoff {xvar_cutoff}\" if xvar_cutoff >= 0 else \"\",\n",
+    "]))\n",
+    "twas_maf_xvar = \" \".join([\n",
+    "    f\"--maf-cutoff {maf_cutoff if maf_cutoff >= 0 else twas_maf_cutoff}\",\n",
+    "    f\"--xvar-cutoff {xvar_cutoff if xvar_cutoff >= 0 else twas_xvar_cutoff}\",\n",
+    "])\n",
     "input: qtl_dataset, for_each = \"region_name\"\n",
     "output: fine_mapping = f\"{cwd:a}/fine_mapping/{name}.{_region_name}.univariate_bvsr.rds\",\n",
     "        twas_weights = f\"{cwd:a}/twas_weights/{name}.{_region_name}.univariate_twas_weights.rds\"\n",
@@ -764,7 +723,7 @@
     "        --methods ${fine_mapping_methods} \\\n",
     "        --coverage ${fine_mapping_coverage} \\\n",
     "        --pip-cutoff-to-skip ${pip_cutoff_to_skip} \\\n",
-    "        ${seed_arg} ${filter_overrides} \\\n",
+    "        ${seed_arg} ${filter_overrides} ${fm_maf_xvar} \\\n",
     "        ${(\"--contexts \" + contexts) if contexts else \"\"} \\\n",
     "        ${(\"--method-args '\" + fine_mapping_method_args + \"'\") if fine_mapping_method_args else \"\"} \\\n",
     "        --output ${_output[\"fine_mapping\"]}\n",
@@ -774,13 +733,11 @@
     "        --gene-id ${_region_name} \\\n",
     "        --cis-window ${cis_window} \\\n",
     "        --methods ${twas_methods} \\\n",
-    "        --min-twas-maf ${min_twas_maf} \\\n",
-    "        --min-twas-xvar ${min_twas_xvar} \\\n",
     "        --max-cv-variants ${max_cv_variants} \\\n",
     "        --cv-folds ${cv_folds} \\\n",
     "        --cv-threads ${cv_threads} \\\n",
     "        --fine-mapping-result ${_output[\"fine_mapping\"]} \\\n",
-    "        ${seed_arg} ${filter_overrides} \\\n",
+    "        ${seed_arg} ${filter_overrides} ${twas_maf_xvar} \\\n",
     "        ${(\"--contexts \" + contexts) if contexts else \"\"} \\\n",
     "        ${(\"--method-args '\" + twas_method_args + \"'\") if twas_method_args else \"\"} \\\n",
     "        --output ${_output[\"twas_weights\"]}"
@@ -994,7 +951,7 @@
     "stop_if(len(regional_data['data']) == 0, f'Either genotype or phenotype data are not available for region {\", \".join(region_name)}.')\n",
     "\n",
     "meta_info = regional_data['meta_info']\n",
-    "input: regional_data[\"data\"], group_by = lambda x: group_by_region(x, regional_data[\"data\"]), group_with = \"meta_info\"\n",
+    "input: regional_data[\"data\"], group_by = lambda x: regional_data[\"data\"], group_with = \"meta_info\"\n",
     "output: f'{cwd:a}/{step_name[:-2]}/{name}.{_meta_info[0]}.mvfsusie_{prior}.rds'\n",
     "task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'\n",
     "R: expand = '${ }', stdout = f\"{_output:n}.stdout\", stderr = f\"{_output:n}.stderr\", container = container, entrypoint = entrypoint\n",
Original file line number	Diff line number	Diff line change
Expand Up		@@ -55,3 +55,5 @@ renovated_code/snakemake/archive/
		renovated_code/snakemake/dryrun/bin/

		dev/

		pixi.lock