diff --git a/code/SoS/data_preprocessing/genotype/genotype_formatting.ipynb b/code/SoS/data_preprocessing/genotype/genotype_formatting.ipynb
index 05b819244..4c87fdbe9 100644
--- a/code/SoS/data_preprocessing/genotype/genotype_formatting.ipynb
+++ b/code/SoS/data_preprocessing/genotype/genotype_formatting.ipynb
@@ -1172,9 +1172,11 @@
     "output: out\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'\n",
+    "    # The paths go in a file: a genome-wide region list is too long for one command-line argument.\n",
+    "    printf '%s\\n' ${_input:r} > \"${_output:n}.data_files.txt\"\n",
     "    Rscript ${modular_script_dir}/data_preprocessing/genotype/genotype_formatting.R \\\n",
     "        --step write_data_list \\\n",
-    "        --data-files \"${\"::\".join([str(x) for x in _input])}\" \\\n",
+    "        --data-files-list \"${_output:n}.data_files.txt\" \\\n",
     "        --ext \"${ext}\" \\\n",
     "        --output \"${_output}\" \\\n",
     "        --numThreads ${numThreads}"
diff --git a/code/SoS/data_preprocessing/phenotype/phenotype_formatting.ipynb b/code/SoS/data_preprocessing/phenotype/phenotype_formatting.ipynb
index 0f773691d..794d00cfc 100644
--- a/code/SoS/data_preprocessing/phenotype/phenotype_formatting.ipynb
+++ b/code/SoS/data_preprocessing/phenotype/phenotype_formatting.ipynb
@@ -611,9 +611,11 @@
     "input: group_by = \"all\"\n",
     "output: f'{cwd}/{name}.{step_name[:-2]}_files.txt'\n",
     "bash: expand= \"${ }\", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', container = container, entrypoint = entrypoint\n",
+    "    # The paths go in a file: a genome-wide region list is too long for the command line.\n",
+    "    printf '%s\\n' ${_input:r} > \"${_output:n}.inputs.txt\"\n",
     "    Rscript ${modular_script_dir}/data_preprocessing/phenotype/phenotype_formatting.R --step phenotype_by_region_2 \\\n",
     "        --cwd \"${cwd}\" \\\n",
-    "        --inputs ${_input} \\\n",
+    "        --inputs-list \"${_output:n}.inputs.txt\" \\\n",
     "        --output \"${_output}\""
    ]
   },
diff --git a/code/SoS/mnm_analysis/mnm_methods/colocboost.ipynb b/code/SoS/mnm_analysis/mnm_methods/colocboost.ipynb
index 2578361c1..23aa15da4 100644
--- a/code/SoS/mnm_analysis/mnm_methods/colocboost.ipynb
+++ b/code/SoS/mnm_analysis/mnm_methods/colocboost.ipynb
@@ -497,7 +497,9 @@
     "parameter: ld_meta_data = path()\n",
     "# Optional region list whose LAST column lists region/gene IDs to analyze.\n",
     "parameter: region_list = path()\n",
-    "# Region/gene IDs to analyze (space separated on the CLI).\n",
+    "# Region/gene IDs to analyze (space separated on the CLI). For thousands of IDs\n",
+    "# use --region-list instead: the IDs reach colocboost_manifest.R as one\n",
+    "# comma-joined argument, which Linux caps at 128 KiB.\n",
     "parameter: region_name = []\n",
     "# Optional per-gene association windows (BED: #chr start end ID) overriding the\n",
     "# default cis_window for the GWAS extraction window.\n",
diff --git a/code/SoS/mnm_analysis/mnm_methods/mnm_regression.ipynb b/code/SoS/mnm_analysis/mnm_methods/mnm_regression.ipynb
index cb2406300..ad779fdd3 100644
--- a/code/SoS/mnm_analysis/mnm_methods/mnm_regression.ipynb
+++ b/code/SoS/mnm_analysis/mnm_methods/mnm_regression.ipynb
@@ -141,7 +141,8 @@
     "- `--customized-association-windows` **`tests/fixtures/qtl_mini/association_windows.bed`**\n",
     "(association windows per region. Without it a fixed window of `--cis-window` is used.)\n",
     "- `--cis-window 1000000` (half-width of the fixed association window used when no window file is given)\n",
-    "- `--region-name ENSG00000283047` (restrict the run to named regions; the union of all regions is used if omitted)\n",
+    "- `--region-name ENSG00000283047` (restrict the run to named regions. Without it or `--region-list`, every region is analyzed: for `susie_twas` and `mnm`, every gene in the phenotype list)\n",
+    "- `--region-list` (a file whose last column names the regions to analyze, for more than fit on the command line; ignored when `--region-name` is given)\n",
     "- `--save-data` (keep the assembled per-region dataset alongside the fine-mapping result)\n",
     "- `--transpose-covariates` (set when the covariate file has samples in columns rather than rows)\n",
     "- `--keep-samples` (a file of sample IDs to retain before fine-mapping)\n",
@@ -665,9 +666,9 @@
     "                        Otherwise, all regions with both genotype and phenotype\n",
     "                        files will be analyzed\n",
     "  --region-name  (as list)\n",
-    "                        Optional: if a region name is provided the analysis\n",
-    "                        would be focused on the union of provides region list\n",
-    "                        and region names\n",
+    "                        Optional: if a region name is provided, the analysis is\n",
+    "                        focused on the named regions and --region-list is\n",
+    "                        ignored.\n",
     "  --keep-samples . (as path)\n",
     "                        Only focus on a subset of samples\n",
     "  --keep-variants . (as path)\n",
@@ -778,7 +779,8 @@
     "      --[no-]drop-indel (default to False)\n",
     "                        Drop indel variants (default keeps them, matching\n",
     "                        QtlDataset()).\n",
-    "  susie_twas:\n",
+    "  susie_twas_1:\n",
+    "  susie_twas_2:\n",
     "    Workflow Options:\n",
     "      --fine-mapping-methods susie\n",
     "      --twas-methods 'susie,mrash,enet,lasso,mcp,scad,l0learn,bayes_r,bayes_c'\n",
@@ -810,7 +812,8 @@
     "                        pipelines, e.g. '{\"susie\":{\"L\":10}}' for fine-mapping or\n",
     "                        '{\"lasso\":{\"nfolds\":10}}' for TWAS.\n",
     "      --twas-method-args ''\n",
-    "  mnm:\n",
+    "  mnm_1:\n",
+    "  mnm_2:\n",
     "    Workflow Options:\n",
     "      --mvsusie-max-iter 200 (as int)\n",
     "      --pip-cutoff-to-skip 0.0 (as float)\n",
@@ -909,8 +912,8 @@
     "# The LAST column of this list will contain the ID of regions to focus on\n",
     "# Otherwise, all regions with both genotype and phenotype files will be analyzed\n",
     "parameter: region_list = path()\n",
-    "# Optional: if a region name is provided \n",
-    "# the analysis would be focused on the union of provides region list and region names\n",
+    "# Optional: if a region name is provided, the analysis is focused on the named\n",
+    "# regions and --region-list is ignored.\n",
     "parameter: region_name = []\n",
     "# Only focus on a subset of samples\n",
     "parameter: keep_samples = path()\n",
@@ -1063,16 +1066,34 @@
    },
    "outputs": [],
    "source": [
-    "[susie_twas]\n",
+    "[susie_twas_1]\n",
+    "# Resolve the gene(s) to analyze -- --region-name, else the last column of\n",
+    "# --region-list, else every gene in the phenotype manifest -- into a manifest via\n",
+    "# region_manifest.R; [susie_twas_2] fans out over its rows. The names go in a\n",
+    "# file (thousands are too long for the command line), and the files\n",
+    "# region_manifest.R reads are inputs, so editing one reruns this step.\n",
+    "input: [phenoFile[0]] + ([region_list] if region_list != path() else [])\n",
+    "output: f\"{cwd:a}/fine_mapping/{name}.susie_twas_region_manifest.tsv\"\n",
+    "task: trunk_workers = 1, trunk_size = 1, walltime = walltime, mem = mem, cores = numThreads, tags = f\"{step_name}_{_output:bn}\"\n",
+    "bash: expand = \"${ }\", stderr = f\"{_output}.stderr\", stdout = f\"{_output}.stdout\", container = container, entrypoint = entrypoint\n",
+    "    printf '%s\\n' ${\" \".join(region_name)} > \"${_output:n}.region_names.txt\"\n",
+    "    Rscript ${modular_script_dir}/pecotmr_integration/region_manifest.R \\\n",
+    "        --pheno-manifest ${phenoFile[0]} \\\n",
+    "        --gene-mode \\\n",
+    "        --cis-window ${cis_window if cis_window > 0 else 1000000} \\\n",
+    "        --region-name-list \"${_output:n}.region_names.txt\" \\\n",
+    "        ${(\"--region-list \" + str(region_list)) if region_list != path() else \"\"} \\\n",
+    "        --output ${_output}\n",
+    "\n",
+    "[susie_twas_2]\n",
     "# Univariate SuSiE fine-mapping + TWAS weights over the pre-built QtlDataset\n",
-    "# (from the qtl_dataset_construct step), one gene per fan-out unit. The gene(s)\n",
-    "# to analyze are passed on the command line via --region-name; the QtlDataset\n",
-    "# RDS already holds the genotype + per-context phenotypes, so this step neither\n",
-    "# re-reads the phenotype manifest nor rebuilds regional_data. Fine-mapping runs\n",
-    "# first (pecotmr_integration/fine_mapping.R -> fineMappingPipeline); the TWAS\n",
-    "# weights pass then reuses its SuSiE fits via --fine-mapping-result\n",
-    "# (pecotmr_integration/twas_weights.R -> twasWeightsPipeline). Replaces the\n",
-    "# legacy regional_data + susie_twas.R path.\n",
+    "# (from the qtl_dataset_construct step), one gene per fan-out unit: a row of the\n",
+    "# [susie_twas_1] manifest. The QtlDataset RDS already holds the genotype +\n",
+    "# per-context phenotypes, so this step doesn't rebuild regional_data.\n",
+    "# Fine-mapping runs first (pecotmr_integration/fine_mapping.R ->\n",
+    "# fineMappingPipeline); the TWAS weights pass then reuses its SuSiE fits via\n",
+    "# --fine-mapping-result (pecotmr_integration/twas_weights.R ->\n",
+    "# twasWeightsPipeline). Replaces the legacy regional_data + susie_twas.R path.\n",
     "#\n",
     "# Both passes get the same --cv-folds: a cross-validated fine-mapping-derived\n",
     "# weight method (susie) needs each fold's OWN fine-mapping fit, so the\n",
@@ -1111,7 +1132,9 @@
     "# '{\"susie\":{\"L\":10}}' for fine-mapping or '{\"lasso\":{\"nfolds\":10}}' for TWAS.\n",
     "parameter: fine_mapping_method_args = \"\"\n",
     "parameter: twas_method_args = \"\"\n",
-    "fail_if(len(region_name) == 0, \"susie_twas: pass --region-name <gene-id> [<gene-id> ...] to choose which gene(s) to analyze.\")\n",
+    "import csv\n",
+    "manifest = path(f\"{cwd:a}/fine_mapping/{name}.susie_twas_region_manifest.tsv\")\n",
+    "jobs = list(csv.DictReader(open(manifest), delimiter='\\t'))\n",
     "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",
@@ -1134,14 +1157,14 @@
     "    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",
+    "input: qtl_dataset, for_each = \"jobs\"\n",
+    "output: fine_mapping = f\"{cwd:a}/fine_mapping/{name}.{_jobs['region_name']}.univariate_bvsr.rds\",\n",
+    "        twas_weights = f\"{cwd:a}/twas_weights/{name}.{_jobs['region_name']}.univariate_twas_weights.rds\"\n",
     "task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f\"{step_name}_{_output[0]:bnn}\"\n",
     "bash: expand = \"${ }\", stderr = f\"{_output[0]:nn}.susie_twas.stderr\", stdout = f\"{_output[0]:nn}.susie_twas.stdout\", container = container, entrypoint = entrypoint\n",
     "    Rscript ${modular_script_dir}/pecotmr_integration/fine_mapping.R \\\n",
     "        --qtl-dataset ${_input} \\\n",
-    "        --gene-id ${_region_name} \\\n",
+    "        --gene-id ${_jobs['region_name']} \\\n",
     "        --cis-window ${cis_window} \\\n",
     "        --methods ${fine_mapping_methods} \\\n",
     "        --coverage ${fine_mapping_coverage} \\\n",
@@ -1155,7 +1178,7 @@
     "\n",
     "    Rscript ${modular_script_dir}/pecotmr_integration/twas_weights.R \\\n",
     "        --qtl-dataset ${_input} \\\n",
-    "        --gene-id ${_region_name} \\\n",
+    "        --gene-id ${_jobs['region_name']} \\\n",
     "        --cis-window ${cis_window} \\\n",
     "        --methods ${twas_methods} \\\n",
     "        --max-cv-variants ${max_cv_variants} \\\n",
@@ -1177,10 +1200,27 @@
    },
    "outputs": [],
    "source": [
-    "[mnm]\n",
+    "[mnm_1]\n",
+    "# Resolve the gene(s) to analyze into a manifest via region_manifest.R, as\n",
+    "# [susie_twas_1] does; [mnm_2] fans out over its rows.\n",
+    "input: [phenoFile[0]] + ([region_list] if region_list != path() else [])\n",
+    "output: f\"{cwd:a}/multivariate_fine_mapping/{name}.mnm_region_manifest.tsv\"\n",
+    "task: trunk_workers = 1, trunk_size = 1, walltime = walltime, mem = mem, cores = numThreads, tags = f\"{step_name}_{_output:bn}\"\n",
+    "bash: expand = \"${ }\", stderr = f\"{_output}.stderr\", stdout = f\"{_output}.stdout\", container = container, entrypoint = entrypoint\n",
+    "    printf '%s\\n' ${\" \".join(region_name)} > \"${_output:n}.region_names.txt\"\n",
+    "    Rscript ${modular_script_dir}/pecotmr_integration/region_manifest.R \\\n",
+    "        --pheno-manifest ${phenoFile[0]} \\\n",
+    "        --gene-mode \\\n",
+    "        --cis-window ${cis_window if cis_window > 0 else 1000000} \\\n",
+    "        --region-name-list \"${_output:n}.region_names.txt\" \\\n",
+    "        ${(\"--region-list \" + str(region_list)) if region_list != path() else \"\"} \\\n",
+    "        --output ${_output}\n",
+    "\n",
+    "[mnm_2]\n",
     "# Multivariate cross-context fine-mapping with mvSuSiE over the pre-built\n",
-    "# QtlDataset (from qtl_dataset_construct), one gene per fan-out unit. Each gene's\n",
-    "# contexts are fit jointly (jointSpecification=\"context\"). The canonical mixture\n",
+    "# QtlDataset (from qtl_dataset_construct), one gene per fan-out unit: a row of\n",
+    "# the [mnm_1] manifest. Each gene's contexts are fit jointly\n",
+    "# (jointSpecification=\"context\"). The canonical mixture\n",
     "# prior is used by default; pass --prior-twas-weights (a TwasWeights RDS from a\n",
     "# preceding `twas_weights.R --methods mrmash` run) for the mr.mash data-driven\n",
     "# reweighted prior. Replaces the legacy mnm.R + multivariate_analysis_pipeline.\n",
@@ -1194,16 +1234,18 @@
     "parameter: prior_weights_min = 1e-10\n",
     "# Optional JSON of extra mvsusie kwargs (overrides the max_iter default), e.g. '{\"mvsusie\":{\"max_iter\":500}}'.\n",
     "parameter: mvsusie_method_args = \"\"\n",
-    "fail_if(len(region_name) == 0, \"mnm: pass --region-name <gene-id> [<gene-id> ...] to choose which gene(s) to analyze.\")\n",
+    "import csv\n",
+    "manifest = path(f\"{cwd:a}/multivariate_fine_mapping/{name}.mnm_region_manifest.tsv\")\n",
+    "jobs = list(csv.DictReader(open(manifest), delimiter='\\t'))\n",
     "qtl_dataset = path(f\"{cwd:a}/qtl_dataset/{name}.qtl_dataset.rds\")\n",
     "mvsusie_args = mvsusie_method_args if mvsusie_method_args else f'{{\"mvsusie\":{{\"max_iter\":{mvsusie_max_iter}}}}}'\n",
-    "input: qtl_dataset, for_each = \"region_name\"\n",
-    "output: f\"{cwd:a}/multivariate_fine_mapping/{name}.{_region_name}.multicontext_bvsr.rds\"\n",
+    "input: qtl_dataset, for_each = \"jobs\"\n",
+    "output: f\"{cwd:a}/multivariate_fine_mapping/{name}.{_jobs['region_name']}.multicontext_bvsr.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",
     "    Rscript ${modular_script_dir}/pecotmr_integration/fine_mapping.R \\\n",
     "        --qtl-dataset ${_input} \\\n",
-    "        --gene-id ${_region_name} \\\n",
+    "        --gene-id ${_jobs['region_name']} \\\n",
     "        --cis-window ${cis_window} \\\n",
     "        --methods mvsusie \\\n",
     "        --joint-specification context \\\n",
@@ -1232,17 +1274,20 @@
     "# Resolve per-locus regions (region_id -> chr:start-end) into a manifest via\n",
     "# region_manifest.R for the cross-gene (jointSpecification=\"trait\") mvSuSiE fit.\n",
     "# Loci come from --customized-association-windows (or --region-list); restrict\n",
-    "# with --region-name. Falls back to per-gene cis windows from the manifest.\n",
-    "input: None\n",
+    "# with --region-name. Falls back to per-gene cis windows from the manifest. The\n",
+    "# names go in a file (thousands are too long for the command line), and the\n",
+    "# files region_manifest.R reads are inputs, so editing one reruns this step.\n",
+    "input: [phenoFile[0]] + [x for x in (customized_association_windows, region_list) if x != path()]\n",
     "output: f\"{cwd:a}/multivariate_fine_mapping/{name}.mnm_genes_region_manifest.tsv\"\n",
     "task: trunk_workers = 1, trunk_size = 1, walltime = walltime, mem = mem, cores = numThreads, tags = f\"{step_name}_{_output:bn}\"\n",
     "bash: expand = \"${ }\", stderr = f\"{_output}.stderr\", stdout = f\"{_output}.stdout\", container = container, entrypoint = entrypoint\n",
+    "    printf '%s\\n' ${\" \".join(region_name)} > \"${_output:n}.region_names.txt\"\n",
     "    Rscript ${modular_script_dir}/pecotmr_integration/region_manifest.R \\\n",
     "        --pheno-manifest ${phenoFile[0]} \\\n",
     "        --cis-window ${cis_window if cis_window > 0 else 1000000} \\\n",
-    "        ${(\"--customized-association-windows \" + str(customized_association_windows)) if customized_association_windows.is_file() else \"\"} \\\n",
-    "        ${(\"--region-name \" + \",\".join(region_name)) if len(region_name) > 0 else \"\"} \\\n",
-    "        ${(\"--region-list \" + str(region_list)) if region_list.is_file() else \"\"} \\\n",
+    "        ${(\"--customized-association-windows \" + str(customized_association_windows)) if customized_association_windows != path() else \"\"} \\\n",
+    "        --region-name-list \"${_output:n}.region_names.txt\" \\\n",
+    "        ${(\"--region-list \" + str(region_list)) if region_list != path() else \"\"} \\\n",
     "        --output ${_output}\n",
     "\n",
     "[mnm_genes_2]\n",
@@ -1296,16 +1341,19 @@
     "# functional SuSiE fit (one locus = one fsusie unit; fsusieR::susiF jointly fits\n",
     "# the multi-trait function within each context). Loci come from\n",
     "# --customized-association-windows (or --region-list); restrict with --region-name.\n",
-    "input: None\n",
+    "# The names go in a file and the files region_manifest.R reads are inputs, as in\n",
+    "# [mnm_genes_1].\n",
+    "input: [phenoFile[0]] + [x for x in (customized_association_windows, region_list) if x != path()]\n",
     "output: f\"{cwd:a}/fsusie/{name}.fsusie_region_manifest.tsv\"\n",
     "task: trunk_workers = 1, trunk_size = 1, walltime = walltime, mem = mem, cores = numThreads, tags = f\"{step_name}_{_output:bn}\"\n",
     "bash: expand = \"${ }\", stderr = f\"{_output}.stderr\", stdout = f\"{_output}.stdout\", container = container, entrypoint = entrypoint\n",
+    "    printf '%s\\n' ${\" \".join(region_name)} > \"${_output:n}.region_names.txt\"\n",
     "    Rscript ${modular_script_dir}/pecotmr_integration/region_manifest.R \\\n",
     "        --pheno-manifest ${phenoFile[0]} \\\n",
     "        --cis-window ${cis_window if cis_window > 0 else 1000000} \\\n",
-    "        ${(\"--customized-association-windows \" + str(customized_association_windows)) if customized_association_windows.is_file() else \"\"} \\\n",
-    "        ${(\"--region-name \" + \",\".join(region_name)) if len(region_name) > 0 else \"\"} \\\n",
-    "        ${(\"--region-list \" + str(region_list)) if region_list.is_file() else \"\"} \\\n",
+    "        ${(\"--customized-association-windows \" + str(customized_association_windows)) if customized_association_windows != path() else \"\"} \\\n",
+    "        --region-name-list \"${_output:n}.region_names.txt\" \\\n",
+    "        ${(\"--region-list \" + str(region_list)) if region_list != path() else \"\"} \\\n",
     "        --output ${_output}\n",
     "\n",
     "[fsusie_2]\n",
@@ -1376,8 +1424,10 @@
     "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",
+    "    # The paths go in a file: a genome-wide list is too long for the command line.\n",
+    "    printf '%s\\n' ${bed_inputs:r} > \"${_output[0]:nn}.input.txt\"\n",
     "    Rscript ${modular_script_dir}/pecotmr_integration/fine_mapping_bed_merge.R \\\n",
-    "        --input ${bed_inputs:r} \\\n",
+    "        --input-list \"${_output[0]:nn}.input.txt\" \\\n",
     "        --output ${_output[0]:r}"
    ],
    "id": "fsusie-export"
@@ -1448,7 +1498,7 @@
     "| Step | Problem | Possible Reason | Solution |\n",
     "| --- | --- | --- | --- |\n",
     "| susie_twas | A region is skipped / no window found | 4th-column `ID` in the association-windows file does not match the phenotype list `ID` | Make the `ID` columns match exactly. |\n",
-    "| susie_twas | `susie_twas: pass --region-name <gene-id>` | No region selected | Pass `--region-name <gene-id>` (e.g. `ENSG00000130538`) or a `--region-list`. |\n",
+    "| susie_twas | `Gene '<id>' not found in --pheno-manifest.` | A `--region-name` / `--region-list` ID is not in the phenotype list | Use IDs from the phenotype list's `ID` column. |\n",
     "| susie_twas | `ERROR: One of the local workers has been killed` | Out-of-memory: too many regions in one process (`-j4`, or all 10 serially) | Use `-j1`; re-run heavy regions one at a time with a single-region phenotype list. |\n",
     "| mnm | `Target unavailable: ...qtl_dataset.rds` | `qtl_dataset_construct` not run first | Chain it: `qtl_dataset_construct+susie_twas`, then `mnm`. |\n",
     "| mnm | `--phenotype-manifest missing required column(s): cond` | Single-context manifest used | Use the multi-context manifest (`protocol_example.pheno_manifest_context.tsv`). |\n",
diff --git a/code/SoS/mnm_analysis/mnm_postprocessing.ipynb b/code/SoS/mnm_analysis/mnm_postprocessing.ipynb
index c8f4bc6b0..b88557817 100644
--- a/code/SoS/mnm_analysis/mnm_postprocessing.ipynb
+++ b/code/SoS/mnm_analysis/mnm_postprocessing.ipynb
@@ -97,7 +97,7 @@
     "\n",
     "The remaining workflows take these additional flags:\n",
     "\n",
-    "- `--tsv_path output_postproc/*.lbf.tsv` (per-region LBF tables to concatenate; `susie_tsv_collapse`)\n",
+    "- `--tsv_path 'output_postproc/*.lbf.tsv'` (per-region LBF tables to concatenate; `susie_tsv_collapse`. Quote the wildcard so SoS expands it: expanded by the shell, a genome-wide set of files is too long for the command line.)\n",
     "- `--annot_tibble <annotation_tibble.tsv>` (annotation tibble used to label variants; `susie_pip_landscape_plot` and `fsusie_extract_effect`)\n",
     "- `--trait_to_select 1` (how many traits the credible-set UpSet diagram shows; `susie_upsetR_cs_plot`)\n",
     "- `--SNP_list <variant_list.txt>` (variants to annotate; `tmp_annotatation_of_snps`)\n",
@@ -264,7 +264,7 @@
    "source": [
     "sos run pipeline/mnm_postprocessing.ipynb susie_tsv_collapse \\\n",
     "    --cwd output_postproc --study <study_name> \\\n",
-    "    --tsv_path output_postproc/*.lbf.tsv \\\n",
+    "    --tsv_path 'output_postproc/*.lbf.tsv' \\\n",
     "    --container '' -j1"
    ]
   },
@@ -276,7 +276,7 @@
    },
    "source": [
     "### Diagnostic plots\n",
-    "`susie_pip_landscape_plot` draws the PIP landscape; `susie_upsetR_plot` and `susie_upsetR_cs_plot` draw UpSet diagrams of credible-set overlap. The landscape plot needs an annotation tibble; `--trait_to_select` picks how many traits the credible-set diagram shows."
+    "`susie_pip_landscape_plot` draws the PIP landscape; `susie_upsetR_plot` and `susie_upsetR_cs_plot` draw UpSet diagrams of credible-set overlap. The landscape plot needs an annotation tibble; `--trait_to_select` picks how many traits the credible-set diagram shows. Quote the `--rds_path` wildcard so SoS expands it: expanded by the shell, a genome-wide set of files is too long for the command line."
    ]
   },
   {
@@ -300,13 +300,13 @@
    "source": [
     "sos run pipeline/mnm_postprocessing.ipynb susie_pip_landscape_plot \\\n",
     "    --cwd output_postproc --study <study_name> \\\n",
-    "    --rds_path <fine_mapping_dir>/*.fine_mapping.rds \\\n",
+    "    --rds_path '<fine_mapping_dir>/*.fine_mapping.rds' \\\n",
     "    --annot_tibble <annotation_tibble.tsv> \\\n",
     "    --container '' -j1\n",
     "\n",
     "sos run pipeline/mnm_postprocessing.ipynb susie_upsetR_cs_plot \\\n",
     "    --cwd output_postproc --study <study_name> \\\n",
-    "    --rds_path <fine_mapping_dir>/*.fine_mapping.rds \\\n",
+    "    --rds_path '<fine_mapping_dir>/*.fine_mapping.rds' \\\n",
     "    --trait_to_select 1 \\\n",
     "    --container '' -j1"
    ]
@@ -1098,8 +1098,10 @@
     "output: f'{cwd}/{name}.UpSetR.pdf'\n",
     "task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'\n",
     "bash: expand = '${ }', stdout = f\"{_output:n}.stdout\", stderr = f\"{_output:n}.stderr\", container = container, entrypoint = entrypoint\n",
+    "    # The paths go in a file: a genome-wide list is too long for the command line.\n",
+    "    printf '%s\\n' ${_input:r} > \"${_output:n}.input.txt\"\n",
     "    Rscript ${modular_script_dir}/pecotmr_integration/fine_mapping_upset.R \\\n",
-    "        --input ${\" \".join([str(x) for x in _input])} \\\n",
+    "        --input-list \"${_output:n}.input.txt\" \\\n",
     "        --output ${_output}\n"
    ]
   },
@@ -1123,8 +1125,10 @@
     "output: f'{cwd}/{name}.UpSetR_cs.pdf'\n",
     "task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'\n",
     "bash: expand = '${ }', stdout = f\"{_output:n}.stdout\", stderr = f\"{_output:n}.stderr\", container = container, entrypoint = entrypoint\n",
+    "    # The paths go in a file: a genome-wide list is too long for the command line.\n",
+    "    printf '%s\\n' ${_input:r} > \"${_output:n}.input.txt\"\n",
     "    Rscript ${modular_script_dir}/pecotmr_integration/fine_mapping_upset.R \\\n",
-    "        --input ${\" \".join([str(x) for x in _input])} \\\n",
+    "        --input-list \"${_output:n}.input.txt\" \\\n",
     "        --output ${_output}\n"
    ]
   },
@@ -1505,8 +1509,10 @@
     "task: trunk_workers = 1, walltime = '1h', trunk_size = 1, mem = '16G', cores = 1, tags = f'{_output:bn}'\n",
     "bash: expand = \"${ }\", container = container, stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout'\n",
     "    if [ -e \"${_output:ad}/${name}_cache/\" ]; then\n",
-    "        sed 's/^chr/#chr/' `ls ${_output:ad}/${name}_cache/*tsv |head -n1` | head -n 1 > ${_output:an}.temp\n",
-    "        tail -n +2 -q ${_output:ad}/${name}_cache/*.tsv >> ${_output:an}.temp\n",
+    "        # find + xargs, not a glob: a genome-wide run has one cache table per region, too many for one command line.\n",
+    "        first_table=$(find \"${_output:ad}/${name}_cache/\" -maxdepth 1 -type f -name '*tsv' | sort | head -n 1)\n",
+    "        sed 's/^chr/#chr/' \"$first_table\" | head -n 1 > ${_output:an}.temp\n",
+    "        find \"${_output:ad}/${name}_cache/\" -maxdepth 1 -type f -name '*.tsv' -print0 | sort -z | xargs -0 -r tail -n +2 -q >> ${_output:an}.temp\n",
     "        error_files=$(find \"${_output:ad}/${name}_cache/\" -type f -name \"*_error\")\n",
     "\n",
     "        if [[ -n $error_files ]]; then\n",
@@ -1625,8 +1631,10 @@
     "parameter: remove_cache = False\n",
     "output: f\"{cache_path:d}/{output_file}\" \n",
     "bash: expand = \"${ }\", container = container, stderr = f'{_input:n}.stderr', stdout = f'{_input:n}.stdout'   \n",
-    "    head -n 1 -q ${cache_path}/*.tsv | sed 's/^chr/#chr/' | head -n 1 > ${_output}\n",
-    "    tail -n +2 -q ${cache_path}/*.tsv >> ${_output}\n",
+    "    # find + xargs, not a glob: one cache table per region is too many for one command line.\n",
+    "    first_table=$(find \"${cache_path}\" -maxdepth 1 -type f -name '*.tsv' | sort | head -n 1)\n",
+    "    head -n 1 \"$first_table\" | sed 's/^chr/#chr/' > ${_output}\n",
+    "    find \"${cache_path}\" -maxdepth 1 -type f -name '*.tsv' -print0 | sort -z | xargs -0 -r tail -n +2 -q >> ${_output}\n",
     "    if ${\"true\" if remove_cache else \"false\"}; then\n",
     "        rm -rf ${cache_path}\n",
     "    fi"
@@ -1708,8 +1716,10 @@
     "output: f\"{cwd}/summary/{name}.{qtl_type}{_tag}.top_loci.bed.gz\"\n",
     "task: trunk_workers = 1, walltime = '2h', trunk_size = 1, mem = '8G', cores = 1, tags = f'{_output:bn}'\n",
     "bash: expand = \"${ }\", container = container, stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', entrypoint=entrypoint\n",
+    "    # The paths go in a file: a genome-wide list is too long for the command line.\n",
+    "    printf '%s\\n' ${_input:r} > \"${_output:nn}.input.txt\"\n",
     "    Rscript ${modular_script_dir}/pecotmr_integration/fine_mapping_bed_merge.R \\\n",
-    "        --input ${_input:r} \\\n",
+    "        --input-list \"${_output:nn}.input.txt\" \\\n",
     "        --output ${_output:r}"
    ]
   },
diff --git a/code/SoS/molecular_phenotypes/QC/splicing_normalization.ipynb b/code/SoS/molecular_phenotypes/QC/splicing_normalization.ipynb
index 3a9ee8cac..8d0309de6 100644
--- a/code/SoS/molecular_phenotypes/QC/splicing_normalization.ipynb
+++ b/code/SoS/molecular_phenotypes/QC/splicing_normalization.ipynb
@@ -554,7 +554,10 @@
     "output: f'{cwd}/{leafcutter_version}/{dataset}_intron_usage_perind.junc'\n",
     "bash: container=container, expand= \"${ }\", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stderr' \n",
     "    cd ${cwd}\n",
-    "    realpath  ${junc_path}/*${file_suffix} > ${_output}.tmp\n",
+    "    # find, not a glob: a cohort's junction files are too many arguments for\n",
+    "    # realpath. Hidden files are skipped, as the glob skipped them.\n",
+    "    find \"${junc_path}/\" -mindepth 1 -maxdepth 1 -name \"*${file_suffix}\" ! -name '.*' -print0 \\\n",
+    "        | sort -z | xargs -0 -r realpath > ${_output}.tmp\n",
     "    grep -v \"all\" ${_output}.tmp > ${_output}\n",
     "    rm ${_output}.tmp"
    ]
diff --git a/code/SoS/molecular_phenotypes/calling/RNA_calling.ipynb b/code/SoS/molecular_phenotypes/calling/RNA_calling.ipynb
index 800c7a31f..61f75d054 100644
--- a/code/SoS/molecular_phenotypes/calling/RNA_calling.ipynb
+++ b/code/SoS/molecular_phenotypes/calling/RNA_calling.ipynb
@@ -954,11 +954,13 @@
     "input: group_by = \"all\"\n",
     "output: f'{cwd}/{sample_list:n}.trimmed.txt'\n",
     "bash: expand = \"${ }\", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout'\n",
+    "    # The paths go in a file: a large cohort's list is too long for the command line.\n",
+    "    printf '%s\\n' ${_input:r} > \"${_output:n}.input.txt\"\n",
     "    Rscript ${modular_script_dir}/molecular_phenotypes/calling/RNA_calling.R \\\n",
     "        --step fastp_manifest \\\n",
     "        --sample-list \"${sample_list}\" \\\n",
     "        --is-paired-end ${is_paired_end} \\\n",
-    "        --input ${' '.join([f'\"{x}\"' for x in _input])} \\\n",
+    "        --input-list \"${_output:n}.input.txt\" \\\n",
     "        --output \"${_output}\"\n"
    ]
   },
@@ -1255,10 +1257,12 @@
     "output: f'{cwd}/{sample_list:bn}.bam_file_list'\n",
     "task: trunk_workers = 1, trunk_size = 1, walltime = walltime, mem = mem, cores = 1\n",
     "bash: expand= \"${ }\", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout'\n",
+    "    # The paths go in a file: a large cohort's list is too long for the command line.\n",
+    "    printf '%s\\n' ${_input:r} > \"${_output}.input.txt\"\n",
     "    Rscript ${modular_script_dir}/molecular_phenotypes/calling/RNA_calling.R \\\n",
     "        --step star_align_3 \\\n",
     "        --output \"${_output}\" \\\n",
-    "        --input ${\" \".join([str(x) for x in _input])} \\\n",
+    "        --input-list \"${_output}.input.txt\" \\\n",
     "        --sample-id ${\" \".join([str(x) for x in sample_id])} \\\n",
     "        --strand ${\" \".join([str(x) for x in strand])} \\\n",
     "        --var-vcf-file \"${varVCFfile}\"\n"
@@ -1360,11 +1364,13 @@
     "        f'{cwd}/{bam_list:bn}.rnaseqc.metrics.tsv'\n",
     "task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads\n",
     "bash: expand= \"${ }\", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout'\n",
+    "    # The paths go in a file: a large cohort's list is too long for the command line.\n",
+    "    printf '%s\\n' ${_input:r} > \"${cwd}/${bam_list:bn}.rnaseqc.input.txt\"\n",
     "    Rscript ${modular_script_dir}/molecular_phenotypes/calling/RNA_calling.R \\\n",
     "        --step rnaseqc_merge \\\n",
     "        --cwd \"${cwd}\" \\\n",
     "        --name \"${bam_list:bn}\" \\\n",
-    "        --input ${\" \".join([str(x) for x in _input])}\n"
+    "        --input-list \"${cwd}/${bam_list:bn}.rnaseqc.input.txt\"\n"
    ]
   },
   {
@@ -1457,11 +1463,13 @@
     "\n",
     "task: trunk_workers = 1, walltime = walltime, mem = mem, cores = numThreads, trunk_size = job_size\n",
     "bash: expand= \"${ }\", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout'\n",
+    "    # The paths go in a file: a large cohort's list is too long for the command line.\n",
+    "    printf '%s\\n' ${_input:r} > \"${cwd}/${bam_list:bn}.rsem.input.txt\"\n",
     "    Rscript ${modular_script_dir}/molecular_phenotypes/calling/RNA_calling.R \\\n",
     "        --step rsem_call_2 \\\n",
     "        --cwd \"${cwd}\" \\\n",
     "        --name \"${bam_list:bn}\" \\\n",
-    "        --input ${\" \".join([str(x) for x in _input])}\n"
+    "        --input-list \"${cwd}/${bam_list:bn}.rsem.input.txt\"\n"
    ]
   },
   {
diff --git a/code/SoS/multivariate_genome/MASH/mash_posterior.ipynb b/code/SoS/multivariate_genome/MASH/mash_posterior.ipynb
index 3630eafdd..3b8f4f74e 100644
--- a/code/SoS/multivariate_genome/MASH/mash_posterior.ipynb
+++ b/code/SoS/multivariate_genome/MASH/mash_posterior.ipynb
@@ -764,8 +764,10 @@
     "input: group_by = \"all\"\n",
     "output: f\"{cwd}/{name}.posterior_sum.csv\"\n",
     "bash: expand = \"${ }\", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', container = container\n",
+    "    # The paths go in a file: a genome-wide list is too long for the command line.\n",
+    "    printf '%s\\n' ${_input:r} > \"${_output:n}.contrast.txt\"\n",
     "    Rscript ${modular_script_dir}/pecotmr_integration/mash_posterior_contrast_summary.R \\\n",
-    "        --contrast ${_input} \\\n",
+    "        --contrast-list \"${_output:n}.contrast.txt\" \\\n",
     "        --cells ${','.join(cells)} \\\n",
     "        --p-cutoff ${p_cut} \\\n",
     "        --output ${_output}"
@@ -886,8 +888,10 @@
     "input: group_by = \"all\"\n",
     "output: f\"{cwd}/{name}.{step_name}.feature_score_sum.csv\"\n",
     "bash: expand = \"${ }\", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', container = container\n",
+    "    # The paths go in a file: a genome-wide list is too long for the command line.\n",
+    "    printf '%s\\n' ${_input:r} > \"${_output:n}.scores.txt\"\n",
     "    Rscript ${modular_script_dir}/pecotmr_integration/mash_feature_score_merge.R \\\n",
-    "        --scores ${_input} \\\n",
+    "        --scores-list \"${_output:n}.scores.txt\" \\\n",
     "        --output ${_output}"
    ]
   },
@@ -983,8 +987,10 @@
     "input: group_by = \"all\"\n",
     "output: f\"{cwd}/{name}.{step_name}.feature_score_sum.csv\"\n",
     "bash: expand = \"${ }\", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', container = container\n",
+    "    # The paths go in a file: a genome-wide list is too long for the command line.\n",
+    "    printf '%s\\n' ${_input:r} > \"${_output:n}.scores.txt\"\n",
     "    Rscript ${modular_script_dir}/pecotmr_integration/mash_feature_score_merge.R \\\n",
-    "        --scores ${_input} \\\n",
+    "        --scores-list \"${_output:n}.scores.txt\" \\\n",
     "        --output ${_output}"
    ]
   }
diff --git a/code/SoS/multivariate_genome/MASH/mash_preprocessing.ipynb b/code/SoS/multivariate_genome/MASH/mash_preprocessing.ipynb
index a83a15c04..869e80ba8 100644
--- a/code/SoS/multivariate_genome/MASH/mash_preprocessing.ipynb
+++ b/code/SoS/multivariate_genome/MASH/mash_preprocessing.ipynb
@@ -538,14 +538,17 @@
     "with open(fine_mapping_meta) as _f:\n",
     "    _rows = list(csv.DictReader(_f, delimiter = '\\t'))\n",
     "_fmr_paths = [r[finemapping_column] for r in _rows]\n",
-    "_region_ids = \",\".join(str(r[region_column]) for r in _rows)\n",
+    "_region_ids = [str(r[region_column]) for r in _rows]\n",
     "input: _fmr_paths\n",
     "output: f\"{cwd}/{name}.mash_input.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\n",
+    "    # The paths and region ids go in files: genome-wide lists are too long for the command line.\n",
+    "    printf '%s\\n' ${_input:r} > \"${_output:n}.objects.txt\"\n",
+    "    printf '%s\\n' ${\" \".join(repr(x) for x in _region_ids)} > \"${_output:n}.region_ids.txt\"\n",
     "    Rscript ${modular_script_dir}/pecotmr_integration/mash_preprocessing.R \\\n",
-    "        --objects ${_input} \\\n",
-    "        --region-ids ${_region_ids} \\\n",
+    "        --objects-list \"${_output:n}.objects.txt\" \\\n",
+    "        --region-ids-list \"${_output:n}.region_ids.txt\" \\\n",
     "        --coverage ${coverage} \\\n",
     "        --n-random ${n_random} \\\n",
     "        --n-null ${n_null} \\\n",
@@ -622,8 +625,10 @@
     "output: f\"{cwd}/{name}.mash_input.rds\"\n",
     "task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{_output:bn}'\n",
     "bash: expand = \"${ }\", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', container = container\n",
+    "    # The paths go in a file: a genome-wide list is too long for the command line.\n",
+    "    printf '%s\\n' ${_input:r} > \"${_output:n}.objects.txt\"\n",
     "    Rscript ${modular_script_dir}/pecotmr_integration/mash_preprocessing.R \\\n",
-    "        --objects ${_input} \\\n",
+    "        --objects-list \"${_output:n}.objects.txt\" \\\n",
     "        --n-random ${n_random} \\\n",
     "        --n-null ${n_null} \\\n",
     "        --exclude-condition \"${','.join(exclude_condition)}\" \\\n",
diff --git a/code/SoS/pecotmr_integration/SuSiE_enloc.ipynb b/code/SoS/pecotmr_integration/SuSiE_enloc.ipynb
index c114e6635..b0a27d25f 100644
--- a/code/SoS/pecotmr_integration/SuSiE_enloc.ipynb
+++ b/code/SoS/pecotmr_integration/SuSiE_enloc.ipynb
@@ -435,6 +435,8 @@
     "parameter: region_list = path()\n",
     "# Optional: if a region name is provided \n",
     "# the analysis would be focused on the union of provides region list and region names\n",
+    "# For thousands of regions use --region-list instead: the names reach\n",
+    "# enloc_manifest.R as one comma-joined argument, which Linux caps at 128 KiB.\n",
     "parameter: region_name = []\n",
     "# It is required to input the name of the analysis\n",
     "parameter: name = f\"{xqtl_meta_data:bnn}.{gwas_meta_data:bnn}\"\n",
@@ -571,9 +573,12 @@
     "    # The xQTL and GWAS sides are already S4 FineMappingResult collections\n",
     "    # (fine_mapping output); qtl_enrichment.R combines its multi-file inputs\n",
     "    # itself, so the legacy legacy_enloc_finemap_convert.R step is gone.\n",
+    "    # The paths go in files: a condition's genome-wide list is too long for the command line.\n",
+    "    printf '%s\\n' ${_input:r} > \"${_output:n}.qtl_fine_mapping.txt\"\n",
+    "    printf '%s\\n' ${\" \".join(repr(x) for x in _jobs[\"gwas_files\"].split(\",\"))} > \"${_output:n}.gwas_fine_mapping.txt\"\n",
     "    Rscript ${modular_script_dir}/pecotmr_integration/qtl_enrichment.R \\\n",
-    "        --qtl-fine-mapping ${\" \".join([str(x) for x in _input])} \\\n",
-    "        --gwas-fine-mapping ${_jobs[\"gwas_files\"].replace(\",\", \" \")} \\\n",
+    "        --qtl-fine-mapping-list \"${_output:n}.qtl_fine_mapping.txt\" \\\n",
+    "        --gwas-fine-mapping-list \"${_output:n}.gwas_fine_mapping.txt\" \\\n",
     "        --ncore ${numThreads} \\\n",
     "        ${seed_arg} \\\n",
     "        --output ${_output}\n"
diff --git a/code/SoS/pecotmr_integration/twas_ctwas.ipynb b/code/SoS/pecotmr_integration/twas_ctwas.ipynb
index 356761060..41f235e4f 100644
--- a/code/SoS/pecotmr_integration/twas_ctwas.ipynb
+++ b/code/SoS/pecotmr_integration/twas_ctwas.ipynb
@@ -688,14 +688,17 @@
     "      --maxSNP 20000 (as int)\n",
     "      --[no-]multi-group (default to True)\n",
     "      --[no-]skip-assembly (default to False)\n",
-    "      --twas-weights  (as list)\n",
+    "      --twas-weights  paths()\n",
+    "\n",
     "                        TwasWeights source (option c): empty -> the upstream\n",
-    "                        twas step's per-gene weights\n",
-    "                        (the per-gene TwasWeights files in the TWAS manifest);\n",
-    "                        override with a list of prebuilt TwasWeights RDS (e.g.\n",
-    "                        the blessed ctwas weights extracted by the blessed S4\n",
-    "                        ctwas weights) when the toy twas weights carry no\n",
-    "                        signal.\n",
+    "                        twas step's per-gene weights (the per-gene TwasWeights\n",
+    "                        files in the TWAS manifest); override with prebuilt\n",
+    "                        TwasWeights RDS paths (e.g. the blessed S4 ctwas\n",
+    "                        weights) when the toy twas weights carry no signal. A\n",
+    "                        quoted wildcard is expanded.\n",
+    "      --twas-weights-list . (as path)\n",
+    "                        Or a file listing those RDS paths, one per line: a\n",
+    "                        genome-wide list is too long for the command line.\n",
     "  ctwas_2:\n",
     "    Workflow Options:\n",
     "      --[no-]run-param-est (default to False)\n",
@@ -785,6 +788,8 @@
     "parameter: regions = path()\n",
     "# Optional: if a region name is provided \n",
     "# the analysis would be focused on the union of provides region list and region names\n",
+    "# For thousands of regions use --regions instead: the names reach twas_manifest.R\n",
+    "# as one comma-joined argument, which Linux caps at 128 KiB.\n",
     "parameter: region_name = []\n",
     "parameter: name = f\"{xqtl_meta_data:bn}.{gwas_meta_data:bn}\"\n",
     "parameter: container = ''\n",
@@ -1077,29 +1082,30 @@
     "        echo \"NOTE: $n_files weight file(s) across $n_genes gene(s) in region $region_id; using the first file per gene.\" >&2\n",
     "    fi\n",
     "    set -- ${\" \".join([tok for i, (g, f) in enumerate(zip(_jobs[\"genes\"].split(\",\"), _input)) if g not in _jobs[\"genes\"].split(\",\")[:i] for tok in (g, str(f))])}\n",
-    "    W=\"\"; FM=\"\"; OUT=\"\"\n",
+    "    # The per-gene paths go in a table, one row per gene: a chromosome-wide\n",
+    "    # region's lists are too long for the command line.\n",
+    "    gene_table=\"$outdir/${name}.$region_id.twas_genes.tsv\"\n",
+    "    printf 'twas_weights\\tfine_mapping_result\\toutput\\n' > \"$gene_table\"\n",
     "    while [ \"$#\" -ge 2 ]; do\n",
     "        gene=\"$1\"\n",
     "        weight=\"$2\"\n",
     "        shift 2\n",
     "        study=$(basename \"$weight\" | cut -d. -f1)\n",
-    "        W=\"$W,$weight\"\n",
-    "        FM=\"$FM,${fine_mapping_dir}/$study/$study.$gene.fine_mapping.s4.rds\"\n",
-    "        OUT=\"$OUT,$outdir/${name}.$region_id.$gene.twas.rds\"\n",
+    "        printf '%s\\t%s\\t%s\\n' \"$weight\" \\\n",
+    "            \"${fine_mapping_dir}/$study/$study.$gene.fine_mapping.s4.rds\" \\\n",
+    "            \"$outdir/${name}.$region_id.$gene.twas.rds\" >> \"$gene_table\"\n",
     "    done\n",
     "    # One R call per region, so the region GwasSumStats (and its LD reference) is read once for every gene.\n",
     "    Rscript ${modular_script_dir}/pecotmr_integration/twas.R \\\n",
-    "        --twas-weights \"$W\" \\\n",
+    "        --gene-table \"$gene_table\" \\\n",
     "        --gwas-sumstats \"$gss\" \\\n",
     "        --rsq-cutoff ${rsq_cutoff} \\\n",
     "        --rsq-pval-cutoff ${\"Inf\" if rsq_pval_cutoff == float(\"inf\") else rsq_pval_cutoff} \\\n",
     "        --rsq-option ${rsq_option} \\\n",
     "        --rsq-pval-option \"${\",\".join(rsq_pval_option)}\" \\\n",
-    "        --fine-mapping-result \"$FM\" \\\n",
     "        --mr-pval-cutoff ${mr_pval_cutoff} \\\n",
     "        --mr-method ${mr_method} \\\n",
-    "        --mr-pip-cutoff ${mr_pip_cutoff} \\\n",
-    "        --output \"$OUT\""
+    "        --mr-pip-cutoff ${mr_pip_cutoff}"
    ]
   },
   {
@@ -1147,16 +1153,27 @@
     "parameter: numThreads = 4\n",
     "parameter: skip_assembly = False\n",
     "# TwasWeights source (option c): empty -> the upstream twas step's per-gene\n",
-    "# weights (the per-gene TwasWeights files in the TWAS manifest); override with a\n",
-    "# list of prebuilt TwasWeights RDS (e.g. the blessed ctwas weights extracted by\n",
-    "# the blessed S4 ctwas weights) when the toy twas weights carry no signal.\n",
-    "parameter: twas_weights = []\n",
+    "# weights (the per-gene TwasWeights files in the TWAS manifest); override with\n",
+    "# prebuilt TwasWeights RDS paths (e.g. the blessed S4 ctwas weights) when the\n",
+    "# toy twas weights carry no signal. A quoted wildcard is expanded.\n",
+    "parameter: twas_weights = paths()\n",
+    "# Or a file listing those RDS paths, one per line: a genome-wide list is too\n",
+    "# long for the command line.\n",
+    "parameter: twas_weights_list = path()\n",
     "skip_if(skip_assembly == True, \"Skip [ctwas_1] assemble.\")\n",
+    "fail_if(len(twas_weights) > 0 and twas_weights_list != path(), \"Pass --twas-weights or --twas-weights-list, not both.\")\n",
     "jobs = twas_manifest_rows\n",
     "gene_chroms = sorted(set(j[\"chrom\"] for j in jobs if j[\"genes\"]))\n",
     "ctwas_chrom = (f\"chr{int(chrom)}\" if str(chrom).isdigit() else str(chrom)) if str(chrom) else (gene_chroms[0] if len(gene_chroms) == 1 else \"\")\n",
     "stop_if(not ctwas_chrom, f\"Specify --chrom: expected one gene-bearing chromosome, found {gene_chroms}.\")\n",
-    "ctwas_weights = list(twas_weights) if twas_weights else list(dict.fromkeys(f for j in jobs if j[\"chrom\"] == ctwas_chrom for f in j[\"weight_files\"].split(',') if f))\n",
+    "# paths() expands wildcards in strings but not in path objects, so only\n",
+    "# --twas-weights globs; the list file's and the manifest's paths are taken as is.\n",
+    "if twas_weights_list != path():\n",
+    "    ctwas_weights = paths([path(x.strip()) for x in open(twas_weights_list) if x.strip()])\n",
+    "elif twas_weights:\n",
+    "    ctwas_weights = paths([str(w) for w in twas_weights])\n",
+    "else:\n",
+    "    ctwas_weights = paths([path(w) for w in dict.fromkeys(f for j in jobs if j[\"chrom\"] == ctwas_chrom for f in j[\"weight_files\"].split(',') if f)])\n",
     "_chrom_jobs = [j for j in jobs if j[\"chrom\"] == ctwas_chrom]\n",
     "ctwas_studies = _chrom_jobs[0][\"gwas_studies\"].split(',') if _chrom_jobs and _chrom_jobs[0][\"gwas_studies\"] else []\n",
     "ctwas_gwas_files = _chrom_jobs[0][\"gwas_files\"] if _chrom_jobs else \"\"\n",
@@ -1218,9 +1235,11 @@
     "    done\n",
     "\n",
     "    # (3) assemble cTWAS inputs (pecotmr places each gene into its home block).\n",
+    "    # The weight paths go in a file: a genome-scale list overflows the 128 KB per-argument limit.\n",
+    "    printf '%s\\n' ${ctwas_weights:r} > \"${_output:n}.weights.txt\"\n",
     "    Rscript ${modular_script_dir}/pecotmr_integration/ctwas_assemble.R \\\n",
     "        --manifest \"$kept\" \\\n",
-    "        --twas-weights \"${\",\".join([str(w) for w in ctwas_weights])}\" \\\n",
+    "        --twas-weights-list \"${_output:n}.weights.txt\" \\\n",
     "        ${f\"--method {ctwas_method} \" if ctwas_method else \"\"}\\\n",
     "        --twas-weight-cutoff ${twas_weight_cutoff} \\\n",
     "        --cs-min-cor ${cs_min_cor} \\\n",
diff --git a/code/script/data_preprocessing/genotype/genotype_formatting.R b/code/script/data_preprocessing/genotype/genotype_formatting.R
index 0b21b09bf..67796b694 100644
--- a/code/script/data_preprocessing/genotype/genotype_formatting.R
+++ b/code/script/data_preprocessing/genotype/genotype_formatting.R
@@ -5,8 +5,9 @@
 #   ld_by_region_plink_1  package plink's `.ld` matrix + `.bim` variant IDs (the plink
 #                         `--r square0 --make-just-bim` call runs in the notebook cell)
 #                         as an RDS (was a numpy .npz; nothing in-repo reads the .npz).
-#   write_data_list       Filter a "::"-joined file list to non-empty entries and
-#                         write a two-column (#id, #path) manifest.
+#   write_data_list       Filter a "::"-joined file list (or a --data-files-list file,
+#                         one path per line) to non-empty entries and write a
+#                         two-column (#id, #path) manifest.
 #   vcf_gz_summary        Print row/column/header/preview stats for VCF(.gz) files.
 #   file_size_summary     Print file path + human-readable size.
 
@@ -49,8 +50,14 @@ do_ld_by_region <- function(argv) {
 
 do_write_data_list <- function(argv) {
   stopifnot(nzchar(argv$output))
-  files <- Filter(nzchar, strsplit(argv$data_files, "::", fixed = TRUE)[[1]])
-  if (!length(files)) stop("--data-files is required")
+  # A genome-wide region list is too long for one command-line argument, so the
+  # notebook passes it as a file.
+  files <- if (nzchar(argv$data_files_list)) {
+    Filter(nzchar, trimws(read_lines_any(argv$data_files_list)))
+  } else {
+    Filter(nzchar, strsplit(argv$data_files, "::", fixed = TRUE)[[1]])
+  }
+  if (!length(files)) stop("--data-files or --data-files-list is required")
   n <- length(strsplit(argv$ext, ".", fixed = TRUE)[[1]]) + 1
   ids <- character(0); keep <- character(0)
   for (f in files) {
@@ -93,6 +100,9 @@ p <- add_argument(p, "--region-end", help = "ld_by_region: region end bp", defau
 p <- add_argument(p, "--float-type", help = "filename precision tag", default = "16")
 p <- add_argument(p, "--output", help = "output path", default = "")
 p <- add_argument(p, "--data-files", help = "\"::\"-joined file list", default = "")
+p <- add_argument(p, "--data-files-list",
+                  help = "write_data_list: file listing the data files, one per line (instead of --data-files)",
+                  default = "")
 p <- add_argument(p, "--ext", help = "write_data_list: file extension", default = "")
 p <- add_argument(p, "--numThreads", help = "plink threads", default = "8")
 argv <- parse_args(p)
diff --git a/code/script/data_preprocessing/phenotype/phenotype_formatting.R b/code/script/data_preprocessing/phenotype/phenotype_formatting.R
index 026420e01..32e9cc212 100644
--- a/code/script/data_preprocessing/phenotype/phenotype_formatting.R
+++ b/code/script/data_preprocessing/phenotype/phenotype_formatting.R
@@ -46,6 +46,9 @@ parser <- add_argument(parser, "--region", nargs = Inf,
                        type = "character")
 parser <- add_argument(parser, "--inputs", nargs = Inf,
                        help = "[*_2] per-shard output paths to collate", type = "character")
+parser <- add_argument(parser, "--inputs-list",
+                       help = "[*_2] File listing the --inputs paths, one per line (instead of --inputs)",
+                       type = "character", default = "")
 parser <- add_argument(parser, "--output", help = "Explicit output path", type = "character",
                        default = "")
 parser <- add_argument(parser, "--output-files", help = "[*_2] chrom/region file manifest",
@@ -84,6 +87,12 @@ parser <- add_argument(parser, "--numThreads", help = "Number of threads",
 argv <- parse_args(parser)
 
 if (is.na(argv$step)) stop("--step is required")
+# A genome-wide list of per-region shards is too long for the command line, so
+# the notebook passes it as a file.
+if (nzchar(argv$inputs_list)) {
+  argv$inputs <- trimws(vroom::vroom_lines(argv$inputs_list))
+  argv$inputs <- argv$inputs[nzchar(argv$inputs)]
+}
 if (!nzchar(argv$name) && !is.na(argv$phenoFile)) {
   argv$name <- sub("\\.gct(\\.gz)?$", "", basename(argv$phenoFile))
 }
diff --git a/code/script/molecular_phenotypes/calling/RNA_calling.R b/code/script/molecular_phenotypes/calling/RNA_calling.R
index a9b87e1d3..c25412274 100644
--- a/code/script/molecular_phenotypes/calling/RNA_calling.R
+++ b/code/script/molecular_phenotypes/calling/RNA_calling.R
@@ -11,6 +11,9 @@
 #   star_align_3        — STAR outputs -> bam_file_list manifest (filename derivation)
 #   rsem_call_2         — per-sample RSEM results -> 7 merged matrices + .cnt metrics
 #
+# --input-list FILE (one path per line) may replace --input; the notebook uses it
+# because a large cohort's merge inputs are too long for the command line.
+#
 # The three merge steps are faithful ports of RNA_calling.py: numeric columns are
 # carried through verbatim (numerically identical to pandas; the float repr is
 # cosmetic), label columns and the manifest are byte-identical.
@@ -40,6 +43,8 @@ p <- add_argument(p, "--name", type = "character", default = "",
                   help = "[rnaseqc_merge/rsem_call_2] output prefix")
 p <- add_argument(p, "--input", type = "character", nargs = Inf,
                   help = "interleaved input files for the selected step")
+p <- add_argument(p, "--input-list", type = "character", default = "",
+                  help = "file listing the --input files in order, one per line (instead of --input)")
 # star_align_3
 p <- add_argument(p, "--sample-id", type = "character", nargs = Inf,
                   help = "[star_align_3] per-sample IDs")
@@ -59,6 +64,10 @@ argv <- parse_args(p)
 
 # nargs=Inf args are NA when not supplied for the selected step
 for (k in c("input", "sample_id", "strand")) argv[[k]] <- argv[[k]][!is.na(argv[[k]])]
+if (nzchar(argv$input_list)) {
+  argv$input <- trimws(vroom::vroom_lines(argv$input_list))
+  argv$input <- argv$input[nzchar(argv$input)]
+}
 
 is_paired   <- argv$is_paired_end
 source_dir  <- argv$input_dir
diff --git a/code/script/pecotmr_integration/ctwas_assemble.R b/code/script/pecotmr_integration/ctwas_assemble.R
index bf3609468..6d2eb85c1 100644
--- a/code/script/pecotmr_integration/ctwas_assemble.R
+++ b/code/script/pecotmr_integration/ctwas_assemble.R
@@ -18,9 +18,11 @@
 #   --manifest             Block-grid manifest TSV with columns:
 #                            region_id          (string, unique)
 #                            gwas_sumstats_rds  (per-block GwasSumStats RDS)
-#   --twas-weights         Comma-separated FLAT per-gene weight RDS
+#   --twas-weights         Comma-separated FLAT per-gene weight RDS paths
 #                          (TwasWeights or QtlFineMappingResult; each carries
 #                          `traitPos` provenance for placement)
+#   --twas-weights-list    File listing those paths, one per line, instead (a
+#                          genome-scale list is too long for the command line)
 #   --fine-mapping-results Optional comma-separated FineMappingResult RDS used
 #                          only as the CS / PIP rescue-filter source
 #                          (NOT the weight source)
@@ -47,6 +49,9 @@ parser <- add_argument(parser, "--manifest",
 parser <- add_argument(parser, "--twas-weights",
                        help = "Comma-separated FLAT per-gene weight RDS (TwasWeights or QtlFineMappingResult)",
                        type = "character")
+parser <- add_argument(parser, "--twas-weights-list",
+                       help = "File listing the --twas-weights RDS paths, one per line (instead of --twas-weights)",
+                       type = "character", default = "")
 parser <- add_argument(parser, "--fine-mapping-results",
                        help = "Optional comma-separated FineMappingResult RDS (CS/PIP filter source)",
                        type = "character", default = "")
@@ -75,8 +80,10 @@ argv <- parse_args(parser)
 
 split_paths <- function(s) {
   if (is.na(s) || !nzchar(s)) return(character(0))
-  trimws(strsplit(s, ",", fixed = TRUE)[[1L]])
+  p <- trimws(strsplit(s, ",", fixed = TRUE)[[1L]])
+  p[nzchar(p)]
 }
+readListFile <- function(f) { x <- trimws(readLines(f)); x[nzchar(x)] }
 
 manifest <- read.table(argv$manifest, header = TRUE, sep = "\t",
                        stringsAsFactors = FALSE, check.names = FALSE,
@@ -128,9 +135,9 @@ gwasSumStats <- combineGwasSumStats(gwasParts)
 
 # FLAT weight source: an unnamed list of per-gene weight objects. assembleCtwasInputs
 # combines them and places each gene into its home block by `traitPos`.
-weightPaths <- split_paths(argv$twas_weights)
+weightPaths <- if (nzchar(argv$twas_weights_list)) readListFile(argv$twas_weights_list) else split_paths(argv$twas_weights)
 if (length(weightPaths) == 0L)
-  stop("--twas-weights lists no weight RDS paths.")
+  stop("--twas-weights / --twas-weights-list lists no weight RDS paths.")
 weightObjs <- lapply(weightPaths, readRDS)
 
 # Optional precomputed TWAS-Z and the CS/PIP-filter FineMappingResult.
diff --git a/code/script/pecotmr_integration/fine_mapping_bed_merge.R b/code/script/pecotmr_integration/fine_mapping_bed_merge.R
index 20462f03b..83e8c0e63 100644
--- a/code/script/pecotmr_integration/fine_mapping_bed_merge.R
+++ b/code/script/pecotmr_integration/fine_mapping_bed_merge.R
@@ -20,6 +20,8 @@
 #
 # Inputs:
 #   --input <TABLE> [...]  Per-region tables (all sharing one header).
+#   --input-list <PATH>    File listing those tables, one per line, instead (a
+#                          genome-wide list is too long for the command line).
 #   --output <PATH>        Merged output; bgzipped when it ends in .gz/.bgz.
 
 suppressPackageStartupMessages({
@@ -33,18 +35,30 @@ suppressPackageStartupMessages({
 parser <- arg_parser("Merge per-region fine-mapping tables into one sorted table")
 parser <- add_argument(parser, "--input", help = "Per-region table paths",
                        type = "character", nargs = Inf)
+parser <- add_argument(parser, "--input-list",
+                       help = "File listing the per-region table paths, one per line (instead of --input)",
+                       type = "character", default = "")
 parser <- add_argument(parser, "--output", help = "Merged output path",
                        type = "character")
 argv <- parse_args(parser)
 
-inputs <- as.character(argv$input)
+readListFile <- function(f) { x <- trimws(readLines(f)); x[nzchar(x)] }
+inputs <- if (nzchar(argv$input_list)) readListFile(argv$input_list) else as.character(argv$input)
 if (length(inputs) == 0L) stop("--input requires at least one table path.")
 
 # col_types = "c" throughout: the packed ";"-joined columns must survive
-# untouched and the table is re-emitted verbatim. vroom reads the shared
-# header once and skips it in the remaining files.
-merged <- vroom(inputs, delim = "\t", col_types = vroom::cols(.default = "c"),
-                progress = FALSE)
+# untouched and the table is re-emitted verbatim. One eager (altrep = FALSE)
+# read per file: a multi-file vroom() call, or lazy columns, keeps every file
+# open, and a genome-wide merge of thousands of tables then hits the
+# open-file limit.
+tables <- lapply(inputs, vroom, delim = "\t", col_types = vroom::cols(.default = "c"),
+                 progress = FALSE, altrep = FALSE)
+header <- names(tables[[1L]])
+differs <- !vapply(tables, function(t) identical(names(t), header), logical(1))
+if (any(differs))
+  stop("All tables must share the header of ", inputs[[1L]], "; these differ: ",
+       paste(inputs[differs], collapse = ", "))
+merged <- bind_rows(tables)
 
 cols <- names(merged)
 # A half-open BED starts #chr/start/end; anything else came from elsewhere and
diff --git a/code/script/pecotmr_integration/fine_mapping_upset.R b/code/script/pecotmr_integration/fine_mapping_upset.R
index f73b75452..8d5c8e6d1 100644
--- a/code/script/pecotmr_integration/fine_mapping_upset.R
+++ b/code/script/pecotmr_integration/fine_mapping_upset.R
@@ -9,6 +9,9 @@
 #
 # Inputs:
 #   --input <RDS> [<RDS> ...]  One or more FineMappingResult RDS paths.
+#   --input-list <PATH>        File listing those RDS paths, one per line,
+#                              instead (a genome-wide list is too long for
+#                              the command line).
 #   --output <PNG>             Output PNG path.
 #   --max-sets                 Cap on the number of sets to render (UpSet
 #                              gets noisy past ~10). Default 20. Sets are
@@ -25,6 +28,9 @@ parser <- arg_parser("UpSet plot of credible-set variant overlap across FineMapp
 parser <- add_argument(parser, "--input",
                        help = "One or more FineMappingResult RDS paths",
                        type = "character", nargs = Inf)
+parser <- add_argument(parser, "--input-list",
+                       help = "File listing the FineMappingResult RDS paths, one per line (instead of --input)",
+                       type = "character", default = "")
 parser <- add_argument(parser, "--output",
                        help = "Output PNG path", type = "character")
 parser <- add_argument(parser, "--max-sets",
@@ -45,7 +51,8 @@ if (!requireNamespace("UpSetR", quietly = TRUE))
   stop("fine_mapping_upset.R requires the UpSetR package: ",
        "install.packages('UpSetR')")
 
-inputs <- as.character(argv$input)
+readListFile <- function(f) { x <- trimws(readLines(f)); x[nzchar(x)] }
+inputs <- if (nzchar(argv$input_list)) readListFile(argv$input_list) else as.character(argv$input)
 if (length(inputs) == 0L)
   stop("--input requires at least one RDS path.")
 
diff --git a/code/script/pecotmr_integration/mash_feature_score_merge.R b/code/script/pecotmr_integration/mash_feature_score_merge.R
index 1570f7dd8..fd5aa2b68 100644
--- a/code/script/pecotmr_integration/mash_feature_score_merge.R
+++ b/code/script/pecotmr_integration/mash_feature_score_merge.R
@@ -8,6 +8,8 @@
 #
 # Inputs:
 #   --scores f1 [f2 ...]  per-chunk feature-score RDS files.
+#   --scores-list <f>     file listing those RDS files, one per line, instead (a
+#                         genome-wide list is too long for the command line).
 #   --output <path>       output table; ".csv" -> CSV, else RDS.
 
 suppressPackageStartupMessages(library(argparser))
@@ -15,10 +17,13 @@ suppressPackageStartupMessages(library(argparser))
 p <- arg_parser("Merge per-chunk MASH feature-score tables")
 p <- add_argument(p, "--scores", type = "character", nargs = Inf,
                   help = "per-chunk feature-score RDS files")
+p <- add_argument(p, "--scores-list", type = "character", default = "",
+                  help = "file listing the --scores RDS files, one per line (instead of --scores)")
 p <- add_argument(p, "--output", type = "character", help = "output CSV or RDS")
 argv <- parse_args(p)
 
-files <- as.character(argv$scores)
+readListFile <- function(f) { x <- trimws(readLines(f)); x[nzchar(x)] }
+files <- if (nzchar(argv$scores_list)) readListFile(argv$scores_list) else as.character(argv$scores)
 if (length(files) == 0L) stop("--scores requires at least one RDS file.")
 
 merged <- do.call(rbind, lapply(files, function(f) as.data.frame(readRDS(f))))
diff --git a/code/script/pecotmr_integration/mash_posterior_contrast_summary.R b/code/script/pecotmr_integration/mash_posterior_contrast_summary.R
index c6753e182..c51e35a86 100644
--- a/code/script/pecotmr_integration/mash_posterior_contrast_summary.R
+++ b/code/script/pecotmr_integration/mash_posterior_contrast_summary.R
@@ -10,6 +10,8 @@
 #
 # Inputs:
 #   --contrast f1 [...]  per-region contrast RDS files (from mash_posterior_contrast.R).
+#   --contrast-list <f>  file listing those RDS files, one per line, instead (a
+#                        genome-wide list is too long for the command line).
 #   --cells c1,c2,...     conditions (columns) whose pairwise contrasts to tally.
 #   --p-cutoff X         significance cutoff. Default 1e-5.
 #   --output <CSV>       output summary CSV (rows: n_sig_snp / n_snp /
@@ -22,6 +24,8 @@ suppressPackageStartupMessages({
 p <- arg_parser("Summarize posterior-contrast significance across regions")
 p <- add_argument(p, "--contrast", type = "character", nargs = Inf,
                   help = "per-region contrast RDS files")
+p <- add_argument(p, "--contrast-list", type = "character", default = "",
+                  help = "file listing the --contrast RDS files, one per line (instead of --contrast)")
 p <- add_argument(p, "--cells", type = "character",
                   help = "comma-separated conditions")
 p <- add_argument(p, "--p-cutoff", type = "numeric", default = 1e-5,
@@ -30,7 +34,8 @@ p <- add_argument(p, "--output", type = "character", help = "output summary CSV"
 argv <- parse_args(p)
 
 cells <- trimws(strsplit(argv$cells, ",", fixed = TRUE)[[1L]])
-files <- as.character(argv$contrast)
+readListFile <- function(f) { x <- trimws(readLines(f)); x[nzchar(x)] }
+files <- if (nzchar(argv$contrast_list)) readListFile(argv$contrast_list) else as.character(argv$contrast)
 if (length(files) == 0L) stop("--contrast requires at least one RDS file.")
 
 pairs <- apply(utils::combn(cells, 2L), 2L, paste, collapse = "_vs_")
diff --git a/code/script/pecotmr_integration/mash_preprocessing.R b/code/script/pecotmr_integration/mash_preprocessing.R
index 49570458c..b4c463923 100644
--- a/code/script/pecotmr_integration/mash_preprocessing.R
+++ b/code/script/pecotmr_integration/mash_preprocessing.R
@@ -20,8 +20,13 @@
 # Inputs:
 #   --objects f1 [f2 ...]   Per-region RDS files (QtlSumStats or
 #                            FineMappingResult). Mixed classes are allowed.
+#   --objects-list <path>   File listing the --objects RDS files, one per line
+#                            (instead of --objects; a genome-wide list is too
+#                            long for the command line).
 #   --region-ids id1,...    Optional comma-separated region labels (one per
 #                            object; default = file basenames).
+#   --region-ids-list <path>  File listing the region labels, one per line
+#                            (instead of --region-ids).
 #   --n-random N            Random rows sampled per object. Default 10.
 #   --n-null N              Null (max|z|<2) rows sampled per object. Default 10.
 #   --exclude-condition c   Comma-separated condition names to drop. Default none.
@@ -45,8 +50,12 @@ suppressPackageStartupMessages({
 p <- arg_parser("Assemble MASH input from per-region S4 objects via mashInput")
 p <- add_argument(p, "--objects", type = "character", nargs = Inf,
                   help = "Per-region RDS files (QtlSumStats or FineMappingResult)")
+p <- add_argument(p, "--objects-list", type = "character", default = "",
+                  help = "File listing the --objects RDS files, one per line (instead of --objects)")
 p <- add_argument(p, "--region-ids", type = "character", default = "",
                   help = "Comma-separated region labels (default: file basenames)")
+p <- add_argument(p, "--region-ids-list", type = "character", default = "",
+                  help = "File listing the region labels, one per line (instead of --region-ids)")
 p <- add_argument(p, "--n-random", type = "integer", default = 10L,
                   help = "Random rows sampled per object")
 p <- add_argument(p, "--n-null", type = "integer", default = 10L,
@@ -82,15 +91,16 @@ readIndependentVariantIds <- function(path) {
   as.character(df[[1L]])
 }
 
-paths <- as.character(argv$objects)
+splitCsv <- function(x) if (nzchar(x)) trimws(strsplit(x, ",", fixed = TRUE)[[1L]]) else character(0)
+readListFile <- function(f) { x <- trimws(readLines(f)); x[nzchar(x)] }
+
+paths <- if (nzchar(argv$objects_list)) readListFile(argv$objects_list) else as.character(argv$objects)
 if (length(paths) == 0L) stop("--objects requires at least one RDS file.")
 missing <- paths[!file.exists(paths)]
 if (length(missing) > 0L)
   stop("--objects file(s) not found: ", paste(missing, collapse = ", "))
 
-splitCsv <- function(x) if (nzchar(x)) trimws(strsplit(x, ",", fixed = TRUE)[[1L]]) else character(0)
-
-regionIds <- splitCsv(argv$region_ids)
+regionIds <- if (nzchar(argv$region_ids_list)) readListFile(argv$region_ids_list) else splitCsv(argv$region_ids)
 if (length(regionIds) == 0L) {
   regionIds <- tools::file_path_sans_ext(basename(paths))
 }
diff --git a/code/script/pecotmr_integration/qtl_enrichment.R b/code/script/pecotmr_integration/qtl_enrichment.R
index 61756259a..6bbe1505c 100644
--- a/code/script/pecotmr_integration/qtl_enrichment.R
+++ b/code/script/pecotmr_integration/qtl_enrichment.R
@@ -12,6 +12,9 @@
 #   --qtl-fine-mapping     Path to S4 QtlFineMappingResult RDS
 #                          (output of fine_mapping.ipynb / fineMappingPipeline)
 #   --gwas-fine-mapping    Path to S4 GwasFineMappingResult RDS
+#   --qtl-fine-mapping-list / --gwas-fine-mapping-list
+#                          File listing those RDS paths, one per line, instead
+#                          (a genome-wide list is too long for the command line)
 #   --num-gwas             Number of GWAS variants per study (pass-through)
 #   --pi-qtl               Optional QTL pi value (pass-through)
 #   --lambda               Pass-through (default 1)
@@ -36,6 +39,12 @@ parser <- add_argument(parser, "--qtl-fine-mapping",
 parser <- add_argument(parser, "--gwas-fine-mapping",
                        help = "Path(s) to S4 GwasFineMappingResult RDS (combined if >1)",
                        type = "character", nargs = Inf)
+parser <- add_argument(parser, "--qtl-fine-mapping-list",
+                       help = "File listing the --qtl-fine-mapping RDS paths, one per line (instead of --qtl-fine-mapping)",
+                       type = "character", default = "")
+parser <- add_argument(parser, "--gwas-fine-mapping-list",
+                       help = "File listing the --gwas-fine-mapping RDS paths, one per line (instead of --gwas-fine-mapping)",
+                       type = "character", default = "")
 parser <- add_argument(parser, "--num-gwas",
                        help = "Number of GWAS variants (per study; pass-through)",
                        type = "numeric", default = NA)
@@ -70,8 +79,11 @@ argv <- parse_args(parser)
   ld <- tryCatch(fmrs[[1L]]@ldReference, error = function(e) NULL)
   do.call(combineFineMappingResults, c(fmrs, list(ldReference = ld)))
 }
-qtlFmr  <- .loadCombine(argv$qtl_fine_mapping)
-gwasFmr <- .loadCombine(argv$gwas_fine_mapping)
+readListFile <- function(f) { x <- trimws(readLines(f)); x[nzchar(x)] }
+qtlFmr  <- .loadCombine(if (nzchar(argv$qtl_fine_mapping_list))
+                          readListFile(argv$qtl_fine_mapping_list) else argv$qtl_fine_mapping)
+gwasFmr <- .loadCombine(if (nzchar(argv$gwas_fine_mapping_list))
+                          readListFile(argv$gwas_fine_mapping_list) else argv$gwas_fine_mapping)
 
 # The impN imputation rounds are drawn in C++ (qtl_enrichment.h), so a
 # main-process set.seed() cannot reach them -- the base seed has to be forwarded,
diff --git a/code/script/pecotmr_integration/region_manifest.R b/code/script/pecotmr_integration/region_manifest.R
index 6b3244411..0f1ba6aa4 100644
--- a/code/script/pecotmr_integration/region_manifest.R
+++ b/code/script/pecotmr_integration/region_manifest.R
@@ -1,10 +1,11 @@
 #!/usr/bin/env Rscript
 # region_manifest.R
 #
-# Resolve mnm_regression's [mnm_genes] / [fsusie] per-locus analysis units into
-# a manifest TSV, so the downstream step can fan out region-mode
-# fineMappingPipeline over its rows with inline csv.DictReader (no notebook-local
-# Python parsing). One row per region (locus).
+# Resolve mnm_regression's analysis units into a manifest TSV, so the downstream
+# step can fan out over its rows with inline csv.DictReader (no notebook-local
+# Python parsing). One row per unit: a locus for [mnm_genes] / [fsusie], which
+# run region-mode fineMappingPipeline over its ld_block, or a gene for
+# [susie_twas] / [mnm], which run gene mode on its region_name.
 #
 # A "region" here is a genomic window that a cross-trait (multi-gene) mvSuSiE or
 # a functional fSuSiE fit runs over; fineMappingPipeline auto-selects the traits
@@ -17,32 +18,49 @@
 #      coordinates +/- --cis-window (single-gene fallback).
 #
 # Inputs:
-#   --pheno-manifest   QtlDataset phenotype manifest TSV (ID, #chr, start, end);
-#                      supplies gene coordinates for the single-gene fallback and
-#                      the default region set when no windows/list are given.
+#   --pheno-manifest   QtlDataset phenotype manifest TSV (ID, #chr, start, end;
+#                      only ID with --gene-mode); supplies gene coordinates for
+#                      the single-gene fallback and the default region set when
+#                      no windows/list are given.
+#   --gene-mode        List genes for a gene-mode fit ([susie_twas] / [mnm]):
+#                      each selected ID must be a --pheno-manifest gene, whose
+#                      coordinates are left empty when the manifest has none;
+#                      --customized-association-windows does not apply.
 #   --customized-association-windows  Optional BED (#chr start end region_id).
-#   --region-name      Optional comma-separated region/gene IDs to restrict to.
-#   --region-list      Optional file whose LAST column lists region IDs.
+#   --region-name      Optional comma-separated region/gene IDs to restrict to;
+#                      --region-list is then ignored.
+#   --region-name-list Optional file of those IDs instead, read like --region-list
+#                      (a genome-wide set is too long for the command line).
+#   --region-list      Optional file whose LAST column lists region IDs; lines
+#                      starting with "#" are comments.
 #   --cis-window       bp added on each side of a gene for the single-gene
 #                      fallback window. Default 1000000.
 #   --output           Output manifest TSV path.
+# A windows or list file that is given must exist.
 #
-# Output TSV columns (one row per region): region_id, chr, start, end, ld_block.
+# Output TSV columns (one row per region): region_id, chr, start, end, ld_block,
+# region_name. region_name is the ID as selected, which gene mode needs verbatim
+# (e.g. a versioned ENSG ID); region_id is that ID made safe for file names.
 
 suppressPackageStartupMessages({
   library(argparser)
 })
 
-parser <- arg_parser("Resolve per-locus regions into a manifest TSV for mnm_genes / fsusie")
+parser <- arg_parser("Resolve mnm_regression's loci or genes into a manifest TSV")
 parser <- add_argument(parser, "--pheno-manifest",
                        help = "QtlDataset phenotype manifest TSV (ID, #chr, start, end)",
                        type = "character")
+parser <- add_argument(parser, "--gene-mode", flag = TRUE,
+                       help = "List genes for a gene-mode fit: IDs must be --pheno-manifest genes, coordinates optional, no windows")
 parser <- add_argument(parser, "--customized-association-windows",
                        help = "Optional BED (#chr start end region_id) of loci",
                        type = "character", default = "")
 parser <- add_argument(parser, "--region-name",
                        help = "Optional comma-separated region/gene IDs to restrict to",
                        type = "character", default = "")
+parser <- add_argument(parser, "--region-name-list",
+                       help = "Optional file of --region-name IDs, read like --region-list (instead of --region-name)",
+                       type = "character", default = "")
 parser <- add_argument(parser, "--region-list",
                        help = "Optional file whose last column lists region IDs",
                        type = "character", default = "")
@@ -66,27 +84,36 @@ id_col    <- intersect(c("ID", "gene_id", "phenotype_id"), names(pm))[1L]
 chr_col   <- intersect(c("#chr", "chrom", "chr"),          names(pm))[1L]
 start_col <- intersect(c("start", "Start"),                names(pm))[1L]
 end_col   <- intersect(c("end", "End"),                    names(pm))[1L]
-if (any(is.na(c(id_col, chr_col, start_col, end_col))))
-  stop("--pheno-manifest needs ID, #chr/chrom, start, end columns (got: ",
-       paste(names(pm), collapse = ", "), ").")
+# Gene mode fits on the gene ID alone, so it needs no coordinates.
+hasCoords <- !anyNA(c(chr_col, start_col, end_col))
+if (is.na(id_col) || !(hasCoords || argv$gene_mode))
+  stop("--pheno-manifest needs ",
+       if (argv$gene_mode) "an ID column" else "ID, #chr/chrom, start, end columns",
+       " (got: ", paste(names(pm), collapse = ", "), ").")
 genes  <- as.character(pm[[id_col]])
-gchr   <- chromAdd(pm[[chr_col]])
-gstart <- suppressWarnings(as.integer(pm[[start_col]]))
-gend   <- suppressWarnings(as.integer(pm[[end_col]]))
 uniqGenes <- unique(genes)
-geneCoord <- lapply(uniqGenes, function(g) {
-  idx <- which(genes == g)
-  list(chr = gchr[idx][[1L]], start = min(gstart[idx], na.rm = TRUE),
-       end = max(gend[idx], na.rm = TRUE))
-})
-names(geneCoord) <- uniqGenes
+geneCoord <- list()
+if (hasCoords) {
+  gchr   <- chromAdd(pm[[chr_col]])
+  gstart <- suppressWarnings(as.integer(pm[[start_col]]))
+  gend   <- suppressWarnings(as.integer(pm[[end_col]]))
+  geneCoord <- lapply(uniqGenes, function(g) {
+    idx <- which(genes == g)
+    list(chr = gchr[idx][[1L]], start = min(gstart[idx], na.rm = TRUE),
+         end = max(gend[idx], na.rm = TRUE))
+  })
+  names(geneCoord) <- uniqGenes
+}
 
 # ----- Customized association windows (loci) --------------------------------
 windows <- list()
 windowOrder <- character(0)
 if (nzchar(argv$customized_association_windows) &&
-    argv$customized_association_windows != "." &&
-    file.exists(argv$customized_association_windows)) {
+    argv$customized_association_windows != ".") {
+  if (argv$gene_mode)
+    stop("--customized-association-windows does not apply with --gene-mode.")
+  if (!file.exists(argv$customized_association_windows))
+    stop("--customized-association-windows not found: ", argv$customized_association_windows)
   caw <- readTableNoHeader(argv$customized_association_windows)
   for (i in seq_len(nrow(caw))) {
     rid <- as.character(caw[[4L]][[i]])
@@ -98,20 +125,19 @@ if (nzchar(argv$customized_association_windows) &&
 }
 
 # ----- Resolve the selected region IDs --------------------------------------
-selected <- character(0)
-rn <- trimws(strsplit(argv$region_name, ",")[[1L]])
-rn <- rn[nzchar(rn)]
+# A list file holds one ID per line, in its last whitespace-separated column.
+readRegionIds <- function(path) {
+  if (!file.exists(path)) stop("Region list not found: ", path)
+  lines <- trimws(vroom::vroom_lines(path))
+  lines <- lines[nzchar(lines) & !startsWith(lines, "#")]
+  unique(vapply(strsplit(lines, "\\s+"), function(parts) parts[[length(parts)]], ""))
+}
+rn <- if (nzchar(argv$region_name_list) && argv$region_name_list != ".")
+  readRegionIds(argv$region_name_list) else splitC(argv$region_name)
 if (length(rn) > 0L) {
   selected <- rn
-} else if (nzchar(argv$region_list) && argv$region_list != "." &&
-           file.exists(argv$region_list)) {
-  for (line in readLines(argv$region_list)) {
-    line <- trimws(line)
-    if (!nzchar(line) || startsWith(line, "#")) next
-    parts <- strsplit(line, "\\s+")[[1L]]
-    selected <- c(selected, parts[[length(parts)]])
-  }
-  selected <- unique(selected)
+} else if (nzchar(argv$region_list) && argv$region_list != ".") {
+  selected <- readRegionIds(argv$region_list)
 } else if (length(windowOrder) > 0L) {
   selected <- windowOrder            # all loci from the windows file
 } else {
@@ -132,6 +158,10 @@ resolve <- function(id) {
     c0 <- geneCoord[[id]]
     return(list(chr = c0$chr, start = max(c0$start - cis, 0L), end = c0$end + cis))
   }
+  if (argv$gene_mode && id %in% uniqGenes)
+    return(list(chr = NA_character_, start = NA_integer_, end = NA_integer_))
+  if (argv$gene_mode)
+    stop("Gene '", id, "' not found in --pheno-manifest.")
   stop("Region/gene '", id, "' not found in --customized-association-windows ",
        "or --pheno-manifest; cannot determine its window.")
 }
@@ -142,7 +172,8 @@ for (id in selected) {
   rows[[length(rows) + 1L]] <- data.frame(
     region_id = gsub("[^A-Za-z0-9_]", "_", id),
     chr       = r$chr, start = r$start, end = r$end,
-    ld_block  = sprintf("%s:%d-%d", r$chr, r$start, r$end),
+    ld_block  = if (is.na(r$chr)) NA_character_ else sprintf("%s:%d-%d", r$chr, r$start, r$end),
+    region_name = id,
     stringsAsFactors = FALSE)
 }
 if (length(rows) == 0L)
diff --git a/code/script/pecotmr_integration/twas.R b/code/script/pecotmr_integration/twas.R
index 44c1f4694..42489f2b7 100644
--- a/code/script/pecotmr_integration/twas.R
+++ b/code/script/pecotmr_integration/twas.R
@@ -12,6 +12,10 @@
 #                          LD reference, is read once for all of them)
 #   --gwas-sumstats        Per-region GwasSumStats RDS (one or more studies)
 #   --fine-mapping-result  Optional per-gene FineMappingResult RDS
+#   --gene-table           TSV with one row per gene and columns twas_weights,
+#                          output and optionally fine_mapping_result, in place of
+#                          the three comma-separated lists (a chromosome-wide
+#                          region's lists are too long for the command line)
 #   --rsq-cutoff           CV-R^2 weight-selection cutoff (legacy rsq_cutoff)
 #   --rsq-pval-cutoff      CV-p-value gate for selection (legacy rsq_pval_cutoff)
 #   --rsq-option           CV r-squared metric: rsq / adj_rsq (legacy rsq_option)
@@ -37,6 +41,9 @@ parser <- add_argument(parser, "--gwas-sumstats",
 parser <- add_argument(parser, "--fine-mapping-result",
                        help = "Optional per-gene FineMappingResult RDS; comma-separated, one per --twas-weights entry",
                        type = "character", default = "")
+parser <- add_argument(parser, "--gene-table",
+                       help = "TSV with one row per gene: twas_weights, output and optionally fine_mapping_result (instead of the three comma-separated lists)",
+                       type = "character", default = "")
 parser <- add_argument(parser, "--mr-pip-cutoff",
                        help = "Pass-through PIP cutoff for MR",
                        type = "numeric", default = 0.5)
@@ -76,14 +83,26 @@ combine_methods <- if (nzchar(argv$combine_methods) && argv$combine_methods != "
 rsq_pval_option <- trimws(strsplit(argv$rsq_pval_option, ",", fixed = TRUE)[[1L]])
 
 # One call may score several genes of the same LD block: --twas-weights,
-# --fine-mapping-result and --output are comma-separated lists, one entry per gene.
+# --fine-mapping-result and --output are comma-separated lists, one entry per
+# gene, or the columns of --gene-table, one row per gene.
 splitList <- function(x) {
   x <- trimws(strsplit(x, ",", fixed = TRUE)[[1L]])
   x[nzchar(x)]
 }
-twPaths <- splitList(argv$twas_weights)
-outPaths <- splitList(argv$output)
-fmPaths <- if (nzchar(argv$fine_mapping_result)) trimws(strsplit(sub("^,", "", argv$fine_mapping_result), ",", fixed = TRUE)[[1L]]) else character(0)
+if (nzchar(argv$gene_table)) {
+  genes <- vroom::vroom(argv$gene_table, delim = "\t", col_types = vroom::cols(.default = "c"),
+                        na = character(), progress = FALSE)
+  missing <- setdiff(c("twas_weights", "output"), names(genes))
+  if (length(missing) > 0L)
+    stop("--gene-table is missing column(s): ", paste(missing, collapse = ", "))
+  twPaths <- genes$twas_weights
+  outPaths <- genes$output
+  fmPaths <- if ("fine_mapping_result" %in% names(genes)) genes$fine_mapping_result else character(0)
+} else {
+  twPaths <- splitList(argv$twas_weights)
+  outPaths <- splitList(argv$output)
+  fmPaths <- if (nzchar(argv$fine_mapping_result)) trimws(strsplit(sub("^,", "", argv$fine_mapping_result), ",", fixed = TRUE)[[1L]]) else character(0)
+}
 if (length(outPaths) != length(twPaths)) stop("--output needs one path per --twas-weights entry.")
 if (length(fmPaths) == 0L) fmPaths <- rep("", length(twPaths))
 if (length(fmPaths) != length(twPaths)) stop("--fine-mapping-result needs one path per --twas-weights entry.")
diff --git a/tests/conftest.py b/tests/conftest.py
index 50931cdd9..5389746a8 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -141,8 +141,8 @@ def ok(script, args, t=400):
 @pytest.fixture(scope="session")
 def ctwas_chain(run_r, repo_root, tmp_path_factory):
     """Run the cTWAS chain once over the blessed 20-block fixture and return the
-    three intermediate RDS ({inputs, est, finemap}); the inputs for the
-    ctwas_assemble / ctwas_est / ctwas_finemap per-script tests."""
+    three intermediate RDS ({inputs, est, finemap}) plus the block-grid manifest;
+    the inputs for the ctwas_assemble / ctwas_est / ctwas_finemap per-script tests."""
     import csv
     r = repo_root
     S = r / "code/script/pecotmr_integration"
@@ -184,4 +184,4 @@ def ok(script, args, t=300):
     ok("ctwas_finemap.R",
        ["--est", est, "--L", "5", "--min-nonsnp-pip", "0.5", "--keep-snps",
         "--ncore", "1", "--output", finemap])
-    return {"inputs": inputs, "est": est, "finemap": finemap}
+    return {"inputs": inputs, "est": est, "finemap": finemap, "manifest": manifest}
diff --git a/tests/fixtures/region_manifest/expected/region_manifest.tsv b/tests/fixtures/region_manifest/expected/region_manifest.tsv
index c1eaff387..f29773dd5 100644
--- a/tests/fixtures/region_manifest/expected/region_manifest.tsv
+++ b/tests/fixtures/region_manifest/expected/region_manifest.tsv
@@ -1,17 +1,17 @@
-region_id	chr	start	end	ld_block
-ENSG00000283047	chr22	10439387	11461338	chr22:10439387-11461338
-ENSG00000130538	chr22	15028191	16029139	chr22:15028191-16029139
-ENSG00000231565	chr22	15111758	16113096	chr22:15111758-16113096
-ENSG00000226474	chr22	15115401	16115578	chr22:15115401-16115578
-ENSG00000235992	chr22	15249155	16250825	chr22:15249155-16250825
-ENSG00000272872	chr22	15323196	16323890	chr22:15323196-16323890
-ENSG00000215270	chr22	15354194	16355243	chr22:15354194-16355243
-ENSG00000229286	chr22	15401790	16401911	chr22:15401790-16401911
-ENSG00000233866	chr22	15414720	16415800	chr22:15414720-16415800
-ENSG00000229658	chr22	15923942	16916492	chr22:15923942-16916492
-ENSG00000235969	chr22	16003303	16916492	chr22:16003303-16916492
-ENSG00000249263	chr22	16153891	16916492	chr22:16153891-16916492
-ENSG00000259271	chr22	16171089	16916492	chr22:16171089-16916492
-ENSG00000237689	chr22	16369477	16916492	chr22:16369477-16916492
-ENSG00000225480	chr22	16397444	16916492	chr22:16397444-16916492
-ENSG00000253691	chr22	16414234	16916492	chr22:16414234-16916492
+region_id	chr	start	end	ld_block	region_name
+ENSG00000283047	chr22	10439387	11461338	chr22:10439387-11461338	ENSG00000283047
+ENSG00000130538	chr22	15028191	16029139	chr22:15028191-16029139	ENSG00000130538
+ENSG00000231565	chr22	15111758	16113096	chr22:15111758-16113096	ENSG00000231565
+ENSG00000226474	chr22	15115401	16115578	chr22:15115401-16115578	ENSG00000226474
+ENSG00000235992	chr22	15249155	16250825	chr22:15249155-16250825	ENSG00000235992
+ENSG00000272872	chr22	15323196	16323890	chr22:15323196-16323890	ENSG00000272872
+ENSG00000215270	chr22	15354194	16355243	chr22:15354194-16355243	ENSG00000215270
+ENSG00000229286	chr22	15401790	16401911	chr22:15401790-16401911	ENSG00000229286
+ENSG00000233866	chr22	15414720	16415800	chr22:15414720-16415800	ENSG00000233866
+ENSG00000229658	chr22	15923942	16916492	chr22:15923942-16916492	ENSG00000229658
+ENSG00000235969	chr22	16003303	16916492	chr22:16003303-16916492	ENSG00000235969
+ENSG00000249263	chr22	16153891	16916492	chr22:16153891-16916492	ENSG00000249263
+ENSG00000259271	chr22	16171089	16916492	chr22:16171089-16916492	ENSG00000259271
+ENSG00000237689	chr22	16369477	16916492	chr22:16369477-16916492	ENSG00000237689
+ENSG00000225480	chr22	16397444	16916492	chr22:16397444-16916492	ENSG00000225480
+ENSG00000253691	chr22	16414234	16916492	chr22:16414234-16916492	ENSG00000253691
diff --git a/tests/notebooks/data_preprocessing/genotype/test_genotype_formatting.py b/tests/notebooks/data_preprocessing/genotype/test_genotype_formatting.py
index 8bfc27d82..d45d64820 100644
--- a/tests/notebooks/data_preprocessing/genotype/test_genotype_formatting.py
+++ b/tests/notebooks/data_preprocessing/genotype/test_genotype_formatting.py
@@ -13,8 +13,11 @@
 import os
 import subprocess
 
+import pytest
+
 from helpers.expected import assert_matches_expected
 
+NB = "pipeline/genotype_formatting.ipynb"
 WORKER = "code/script/data_preprocessing/genotype/genotype_formatting.R"
 SH = "code/script/data_preprocessing/genotype/genotype_formatting.sh"
 EXP = "tests/fixtures/genotype_formatting/expected"
@@ -51,13 +54,20 @@ def _rscript(expr):
     return r.stdout.strip()
 
 
-def test_write_data_list(run_r, repo_root, tmp_path):
+@pytest.mark.parametrize("via_list", [False, True], ids=["inline", "list_file"])
+def test_write_data_list(via_list, run_r, repo_root, tmp_path):
     f1 = tmp_path / "s.chr1.bed"; f1.write_text("x")
     f2 = tmp_path / "s.chr2.bed"; f2.write_text("y")
     f3 = tmp_path / "s.chr3.bed"; f3.write_text("")            # empty -> excluded
     out = tmp_path / "list.tsv"
-    files = "::".join(str(x) for x in (f1, f2, f3))
-    p = run_r(repo_root / WORKER, ["--step", "write_data_list", "--data-files", files,
+    if via_list:
+        # the notebook's form: a genome-wide list is too long for one argument
+        listing = tmp_path / "data_files.txt"
+        listing.write_text("".join(f"{x}\n" for x in (f1, f2, f3)))
+        files = ["--data-files-list", listing]
+    else:
+        files = ["--data-files", "::".join(str(x) for x in (f1, f2, f3))]
+    p = run_r(repo_root / WORKER, ["--step", "write_data_list", *files,
                                    "--ext", "bed", "--output", out])
     assert p.returncode == 0, p.stdout + p.stderr
     lines = out.read_text().splitlines()
@@ -67,6 +77,30 @@ def test_write_data_list(run_r, repo_root, tmp_path):
     assert len(lines) == 3                                     # empty chr3 dropped
 
 
+def test_genotype_by_region_manifest(run_sos, repo_root, tmp_path):
+    """genotype_by_region end to end: one PLINK fileset per region, then
+    write_data_list gathers them into the manifest through a list file (a
+    genome-wide region list is too long for one command-line argument)."""
+    regions = tmp_path / "regions.txt"
+    regions.write_text("22\t16000000\t17000000\tregA\n22\t17000000\t18000000\tregB\n")
+    out = tmp_path / "out"
+    p = run_sos(repo_root / NB, "genotype_by_region",
+                {"genoFile": repo_root / f"{BED}.bed", "region_list": regions,
+                 "window": 0, "cwd": out, "name": "test", "numThreads": 1,
+                 "modular_script_dir": repo_root / "code/script"},
+                cwd=repo_root, timeout=600)
+    assert p.returncode == 0, p.stdout + p.stderr
+    by_region = out / "regions_genotype_by_region"
+    lines = (by_region / "protocol_example.ld_genotype.chr22.genotype_by_region_files.txt"
+             ).read_text().splitlines()
+    assert lines[0] == "#id\t#path"
+    rows = [line.split("\t") for line in lines[1:]]
+    assert [r[0] for r in rows] == ["regA", "regB"]
+    assert [os.path.basename(r[1]) for r in rows] == [
+        "protocol_example.ld_genotype.chr22.regA.bed",
+        "protocol_example.ld_genotype.chr22.regB.bed"]
+
+
 def test_ld_by_region(run_r, repo_root, tmp_path):
     out = tmp_path / "region.float16.rds"
     out_prefix = tmp_path / "region"          # the R step strips .float16.rds -> "region"
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 f99f6a050..24df5c5e8 100644
--- a/tests/notebooks/mnm_analysis/mnm_methods/test_mnm_regression.py
+++ b/tests/notebooks/mnm_analysis/mnm_methods/test_mnm_regression.py
@@ -3,7 +3,9 @@
 This is the website-published entrypoint for the QTL fine-mapping / TWAS-weight
 chain. Running the `qtl_dataset_construct+susie_twas` step exercises the SoS cell
 orchestration end-to-end (param wiring, path resolution, step chaining) that the
-direct script tests can't reach, and asserts the same output S4 shapes.
+direct script tests can't reach, and asserts the same output S4 shapes. The
+gene manifest steps ([susie_twas_1] / [mnm_1]) also run alone, on --region-list
+and with no selection.
 """
 from __future__ import annotations
 
@@ -92,6 +94,65 @@ def test_mnm(run_sos, read_rds, repo_root, qtl_mini, tmp_path):
                             rtol=1e-6, atol=1e-8)
 
 
+def _gene_manifest_params(repo_root, qtl_mini, cwd):
+    """The [global] inputs of the gene manifest steps; no QtlDataset is needed,
+    since only the manifest step runs."""
+    return {
+        "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",
+        "modular_script_dir": repo_root / "code/script",
+    }
+
+
+@pytest.mark.parametrize("workflow, folder", [("susie_twas", "fine_mapping"),
+                                              ("mnm", "multivariate_fine_mapping")])
+def test_gene_manifest_from_region_list(run_sos, repo_root, qtl_mini, tmp_path,
+                                        workflow, folder):
+    """The manifest step takes the genes from --region-list's last column, and
+    reruns when the list is edited: the list is a step input, so SoS sees the
+    change instead of reusing the manifest it saved."""
+    cwd = tmp_path / "mnm"
+    genes = tmp_path / "genes.bed"
+    params = {**_gene_manifest_params(repo_root, qtl_mini, cwd), "region-list": genes}
+    manifest = cwd / folder / f"test_study.{workflow}_region_manifest.tsv"
+    for listed in (["ENSG00000130538", GENE], [GENE]):
+        genes.write_text("#chr\tstart\tend\tID\n"
+                         + "".join(f"chr22\t0\t1\t{gene}\n" for gene in listed))
+        p = run_sos(repo_root / "pipeline/mnm_regression.ipynb", f"{workflow}:1",
+                    params, cwd=repo_root, timeout=300)
+        assert p.returncode == 0, p.stdout + p.stderr
+        rows = [line.split("\t") for line in manifest.read_text().splitlines()]
+        column = rows[0].index("region_name")
+        assert [row[column] for row in rows[1:]] == listed
+
+
+@pytest.mark.parametrize("workflow, folder", [("susie_twas", "fine_mapping"),
+                                              ("mnm", "multivariate_fine_mapping")])
+def test_gene_manifest_defaults_to_every_gene(run_sos, repo_root, qtl_mini, tmp_path,
+                                              workflow, folder):
+    """With neither --region-name nor --region-list, every gene in the phenotype
+    manifest is analyzed. The manifest here has no #chr/start/end columns:
+    qtl_dataset_construct does not need them, and gene mode does not either."""
+    cwd = tmp_path / "mnm"
+    pheno = [line.split("\t") for line in
+             (qtl_mini / "protocol_example.pheno_manifest_context.tsv").read_text().splitlines()]
+    keep = [i for i, column in enumerate(pheno[0]) if column not in ("#chr", "start", "end")]
+    no_coordinates = tmp_path / "pheno_manifest.tsv"
+    no_coordinates.write_text("".join("\t".join(row[i] for i in keep) + "\n" for row in pheno))
+    p = run_sos(repo_root / "pipeline/mnm_regression.ipynb", f"{workflow}:1",
+                {**_gene_manifest_params(repo_root, qtl_mini, cwd), "phenoFile": no_coordinates},
+                cwd=repo_root, timeout=300)
+    assert p.returncode == 0, p.stdout + p.stderr
+    genes = list(dict.fromkeys(row[pheno[0].index("ID")] for row in pheno[1:]))
+    manifest = cwd / folder / f"test_study.{workflow}_region_manifest.tsv"
+    rows = [line.split("\t") for line in manifest.read_text().splitlines()]
+    column = rows[0].index("region_name")
+    assert [row[column] for row in rows[1:]] == genes
+
+
 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
diff --git a/tests/notebooks/mnm_analysis/test_mnm_postprocessing.py b/tests/notebooks/mnm_analysis/test_mnm_postprocessing.py
index 298c8b48e..2ddb356bf 100644
--- a/tests/notebooks/mnm_analysis/test_mnm_postprocessing.py
+++ b/tests/notebooks/mnm_analysis/test_mnm_postprocessing.py
@@ -57,9 +57,29 @@ def test_cis_results_export_and_top_loci(run_sos, read_rds, repo_root, tmp_path)
     assert (cwd / "summary" / "protocol_example.eQTL.top_loci.bed.gz").exists(), p.stdout
 
 
+def test_combine_export_meta(run_sos, repo_root, tmp_path):
+    # Concatenate the per-region cache tables under one header, in file-name order.
+    # The step lists them with find rather than a glob, which overflows the command
+    # line for a genome-wide cache. cwd=tmp_path: the step logs to <cwd>/.stderr.
+    cache = tmp_path / "cache"
+    cache.mkdir()
+    for region, start in (("r2", 200), ("r1", 100)):
+        (cache / f"protocol_example_{region}.tsv").write_text(
+            f"chr\tstart\tend\tregion_id\nchr22\t{start}\t{start + 50}\t{region}\n")
+    p = run_sos(repo_root / NB, "combine_export_meta",
+                dict(study="protocol_example", cache_path=cache, output_file="meta.tsv"),
+                cwd=tmp_path, timeout=300)
+    assert p.returncode == 0, p.stdout + p.stderr
+    assert (tmp_path / "meta.tsv").read_text() == (
+        "#chr\tstart\tend\tregion_id\n"
+        "chr22\t100\t150\tr1\n"
+        "chr22\t200\t250\tr2\n")
+
+
 @pytest.mark.parametrize("wf,pattern", [
     ("susie_pip_landscape_plot", "*.pip_landscape_plot.pdf"),
     ("susie_upsetR_plot", "*.UpSetR.pdf"),
+    ("susie_upsetR_cs_plot", "*.UpSetR_cs.pdf"),
 ])
 def test_plot_wrappers(run_sos, repo_root, tmp_path, wf, pattern):
     fx = repo_root / "tests" / "fixtures" / "mnm_postprocessing"
diff --git a/tests/notebooks/molecular_phenotypes/QC/test_splicing_normalization.py b/tests/notebooks/molecular_phenotypes/QC/test_splicing_normalization.py
index 1ba8abea6..bf12c0901 100644
--- a/tests/notebooks/molecular_phenotypes/QC/test_splicing_normalization.py
+++ b/tests/notebooks/molecular_phenotypes/QC/test_splicing_normalization.py
@@ -1,7 +1,8 @@
 """Notebook tier: splicing_normalization.ipynb workers (leafcutter).
 
 splicing_normalization.R ports the leafcutter prepare_phenotype_table.py logic.
-Fixtures + references are the 20-sample MWE leafcutter outputs.
+Fixtures + references are the 20-sample MWE leafcutter outputs. The [Junc_list]
+step runs through the notebook on the splicing_calling junction fixtures.
 Float phenotype values are compared NUMERICALLY (identical doubles): Python's shortest-repr
 and R's differ cosmetically on ~4% of values in the 16th/17th digit, but the doubles — what
 downstream tools parse — are identical.
@@ -83,3 +84,27 @@ def test_jointcall_samples(run_r, repo_root, tmp_path):
     assert p.returncode == 0, p.stdout + p.stderr
     ids = [ln.split("\t")[0] for ln in out.read_text().splitlines()[1:] if ln.strip()]
     assert set(ids) == {"SAMPLE_001", "SAMPLE_002"}          # .final stripped, SAMPLE_999 dropped
+
+
+def test_junc_list(run_sos, repo_root, tmp_path):
+    # [Junc_list] lists a directory's junction files with find (a cohort's files
+    # are too many arguments for realpath): sorted absolute paths, leaving out
+    # other suffixes and hidden files. The fixtures are copied, not linked:
+    # realpath would resolve a link back to the fixture directory, and the
+    # step's grep -v "all" drops every path containing "all" (splicing_calling).
+    fix = repo_root / "tests/fixtures/splicing_calling/leafcutter"
+    juncs = tmp_path / "juncs"
+    juncs.mkdir()
+    for s in ("SAMPLE_002", "SAMPLE_001"):
+        (juncs / f"{s}.junc.gz").write_bytes((fix / f"{s}.junc.gz").read_bytes())
+    (juncs / ".SAMPLE_003.junc.gz").write_bytes(b"")
+    (juncs / "protocol_example.junction_files.txt").write_text("")
+    cwd = tmp_path / "out"
+    p = run_sos(repo_root / "pipeline/splicing_normalization.ipynb", "Junc_list",
+                {"cwd": cwd, "junc_path": juncs, "file_suffix": "junc.gz",
+                 "dataset": "protocol_example"},
+                cwd=repo_root, timeout=300)
+    assert p.returncode == 0, p.stdout + p.stderr
+    listed = (cwd / "leafcutter2" / "protocol_example_intron_usage_perind.junc").read_text()
+    assert listed.splitlines() == [str((juncs / f"{s}.junc.gz").resolve())
+                                   for s in ("SAMPLE_001", "SAMPLE_002")]
diff --git a/tests/notebooks/molecular_phenotypes/calling/test_rna_calling.py b/tests/notebooks/molecular_phenotypes/calling/test_rna_calling.py
index 53ce34434..a8a3b3516 100644
--- a/tests/notebooks/molecular_phenotypes/calling/test_rna_calling.py
+++ b/tests/notebooks/molecular_phenotypes/calling/test_rna_calling.py
@@ -13,17 +13,47 @@
     match except the machine-dependent File path column.
 Fixtures are the MWE per-sample RSEM results (downsampled to 500 rows) + synthetic
 RNA-SeQC GCTs.
+
+The fastp_trim_adaptor and rnaseqc_call workflows also run through the notebook,
+on two samples that share SAMPLE_001's reads, so the merge steps' hand-off of
+their inputs to RNA_calling.R is covered.
 """
 from __future__ import annotations
 
 import gzip
 import subprocess
+from pathlib import Path
+
+import pytest
 
 WORKER = "code/script/molecular_phenotypes/calling/RNA_calling.R"
 SH = "code/script/molecular_phenotypes/calling/RNA_calling.sh"
 FIX = "tests/fixtures/rna_calling"
 
 
+def _input_args(paths, via_list, tmp_path):
+    """--input, or --input-list: the notebook's form (a large cohort's list is
+    too long for the command line)."""
+    if not via_list:
+        return ["--input", *paths]
+    listing = tmp_path / "inputs.txt"
+    listing.write_text("".join(f"{x}\n" for x in paths))
+    return ["--input-list", listing]
+
+
+def _two_sample_sheet(fix, tmp_path):
+    """The fixture sample sheet, with SAMPLE_002's reads linked to SAMPLE_001's
+    (only SAMPLE_001's are committed)."""
+    data = tmp_path / "data"
+    data.mkdir()
+    for s in ("SAMPLE_001", "SAMPLE_002"):
+        for r in ("R1", "R2"):
+            (data / f"{s}_{r}.fastq.gz").symlink_to(fix / "fastq" / f"SAMPLE_001_{r}.fastq.gz")
+    sheet = data / "protocol_example.rnaseq.fastq.list.txt"
+    sheet.write_text((fix / "protocol_example.rnaseq.fastq.list.txt").read_text())
+    return sheet
+
+
 def _rows(text, sep="\t"):
     return [ln.split(sep) for ln in text.splitlines() if ln.strip()]
 
@@ -48,12 +78,14 @@ def _assert_num_match(got_text, exp_text, label_cols, sep="\t", tol=1e-9):
                 assert abs(float(x) - float(y)) < tol, f"num diff row {i} col {j}: {x} vs {y}"
 
 
-def test_rsem_call_2(run_r, repo_root, tmp_path):
+@pytest.mark.parametrize("via_list", [False, True], ids=["inline", "list_file"])
+def test_rsem_call_2(via_list, run_r, repo_root, tmp_path):
     fix = repo_root / FIX
     inp = []
     for s in ("SAMPLE_001", "SAMPLE_002"):
         inp += [fix / f"{s}.rsem.isoforms.results", fix / f"{s}.rsem.genes.results", fix / f"{s}.rsem.cnt"]
-    p = run_r(repo_root / WORKER, ["--step", "rsem_call_2", "--cwd", tmp_path, "--name", "rsem", "--input", *inp])
+    p = run_r(repo_root / WORKER, ["--step", "rsem_call_2", "--cwd", tmp_path, "--name", "rsem",
+                                   *_input_args(inp, via_list, tmp_path)])
     assert p.returncode == 0, p.stdout + p.stderr
     for f in ("rsem_transcripts_expected_count", "rsem_transcripts_tpm", "rsem_transcripts_fpkm",
               "rsem_transcripts_isopct", "rsem_genes_expected_count", "rsem_genes_tpm", "rsem_genes_fpkm"):
@@ -70,7 +102,8 @@ def test_rsem_call_2(run_r, repo_root, tmp_path):
         assert [float(c) for c in keep(rg)[1:]] == [float(c) for c in keep(re_)[1:]]
 
 
-def test_star_align_3(run_r, repo_root, tmp_path):
+@pytest.mark.parametrize("via_list", [False, True], ids=["inline", "list_file"])
+def test_star_align_3(via_list, run_r, repo_root, tmp_path):
     fix = repo_root / FIX
     inp = []
     for s in ("SAMPLE_001", "SAMPLE_002"):
@@ -78,7 +111,8 @@ def test_star_align_3(run_r, repo_root, tmp_path):
         bw = f"{s}.Aligned.sortedByCoord.out_nowasp_noqc.md.bw"
         inp += ["d0", "d1", coord, "d3", bw, "d5"]
     out = tmp_path / "bam_file_list"
-    p = run_r(repo_root / WORKER, ["--step", "star_align_3", "--output", out, "--input", *inp,
+    p = run_r(repo_root / WORKER, ["--step", "star_align_3", "--output", out,
+                                   *_input_args(inp, via_list, tmp_path),
                                    "--sample-id", "SAMPLE_001", "SAMPLE_002", "--strand", "rf", "rf",
                                    "--var-vcf-file", ""])
     assert p.returncode == 0, p.stdout + p.stderr
@@ -92,12 +126,14 @@ def test_detect_strand(run_r, repo_root):
     assert p.stdout.strip() == (fix / "expected" / "SAMPLE_001.strand.txt").read_text().strip()
 
 
-def test_fastp_manifest(run_r, repo_root, tmp_path):
+@pytest.mark.parametrize("via_list", [False, True], ids=["inline", "list_file"])
+def test_fastp_manifest(via_list, run_r, repo_root, tmp_path):
     fix = repo_root / FIX
     trimmed = [f"out/{s}_{r}.fastq.trimmed.fq.gz" for s in ("SAMPLE_001", "SAMPLE_002") for r in ("R1", "R2")]
     out = tmp_path / "fastq.list.trimmed.txt"
     p = run_r(repo_root / WORKER, ["--step", "fastp_manifest", "--sample-list", fix / "protocol_example.rnaseq.fastq.list.txt",
-                                   "--is-paired-end", 1, "--output", out, "--input", *trimmed])
+                                   "--is-paired-end", 1, "--output", out,
+                                   *_input_args(trimmed, via_list, tmp_path)])
     assert p.returncode == 0, p.stdout + p.stderr
     assert out.read_text() == (fix / "expected" / "fastq.list.trimmed.txt").read_text()
 
@@ -317,16 +353,65 @@ def test_rnaseqc_call(run_sh, repo_root, tmp_path):
     assert (tmp_path / "SAMPLE_001.rnaseqc.metrics.tsv").exists()
 
 
-def test_rnaseqc_merge(run_r, repo_root, tmp_path):
+@pytest.mark.parametrize("via_list", [False, True], ids=["inline", "list_file"])
+def test_rnaseqc_merge(via_list, run_r, repo_root, tmp_path):
     fix = repo_root / FIX
     inp = []
     for s in ("SAMPLE_001", "SAMPLE_002"):
         inp += [fix / f"{s}.rnaseqc.gene_tpm.gct.gz", fix / f"{s}.rnaseqc.gene_reads.gct.gz",
                 fix / f"{s}.rnaseqc.exon_reads.gct.gz", fix / f"{s}.rnaseqc.metrics.tsv"]
-    p = run_r(repo_root / WORKER, ["--step", "rnaseqc_merge", "--cwd", tmp_path, "--name", "rnaseqc", "--input", *inp])
+    p = run_r(repo_root / WORKER, ["--step", "rnaseqc_merge", "--cwd", tmp_path, "--name", "rnaseqc",
+                                   *_input_args(inp, via_list, tmp_path)])
     assert p.returncode == 0, p.stdout + p.stderr
     for f in ("gene_tpm.gct.gz", "gene_readsCount.gct.gz", "exon_readsCount.gct.gz"):
         assert gzip.open(tmp_path / f"rnaseqc.rnaseqc.{f}", "rt").read() == \
             gzip.open(fix / "expected" / f"rnaseqc.rnaseqc.{f}", "rt").read()
     assert (tmp_path / "rnaseqc.rnaseqc.metrics.tsv").read_text() == \
         (fix / "expected" / "rnaseqc.rnaseqc.metrics.tsv").read_text()
+
+
+def test_fastp_trim_adaptor_workflow(run_sos, repo_root, tmp_path):
+    # fastp per sample, then [fastp_trim_adaptor_2] points the sample sheet at the
+    # trimmed reads, whose paths it hands to RNA_calling.R in a list file.
+    fix = repo_root / FIX
+    cwd = tmp_path / "out"
+    p = run_sos(repo_root / "pipeline/RNA_calling.ipynb", "fastp_trim_adaptor",
+                {"sample_list": _two_sample_sheet(fix, tmp_path), "cwd": cwd, "numThreads": 1,
+                 "modular_script_dir": repo_root / "code/script"},
+                cwd=repo_root, timeout=600)
+    assert p.returncode == 0, p.stdout + p.stderr
+    # written under cwd at the sheet's absolute path
+    rows = _rows(next(cwd.rglob("*.trimmed.txt")).read_text())
+    assert rows[0] == ["ID", "fq1", "fq2", "strand", "read_length"]
+    assert [(r[0], Path(r[1]).name, Path(r[2]).name) for r in rows[1:]] == \
+        [(s, f"{s}_R1.fastq.trimmed.fq.gz", f"{s}_R2.fastq.trimmed.fq.gz")
+         for s in ("SAMPLE_001", "SAMPLE_002")]
+
+
+def test_rnaseqc_call_workflow(run_sos, repo_root, tmp_path):
+    # RNA-SeQC per sample, then [rnaseqc_call_2] merges the per-sample tables,
+    # whose paths it hands to RNA_calling.R in a list file. SAMPLE_002's BAM is
+    # linked to SAMPLE_001's.
+    fix = repo_root / FIX
+    bams = tmp_path / "bam"
+    bams.mkdir()
+    bam_list = ["sample_id\tstrand\tcoord_bam_list\tBW_list\tSJ_list\ttrans_bam_list"]
+    for s in ("SAMPLE_001", "SAMPLE_002"):
+        for ext in (".bam", ".bam.bai"):
+            (bams / f"{s}.cord{ext}").symlink_to(fix / "bam" / f"SAMPLE_001.cord{ext}")
+        bam_list.append(f"{s}\trf\t{bams / s}.cord.bam\t-\t-\t-")
+    (tmp_path / "protocol_example.bam_file_list").write_text("\n".join(bam_list) + "\n")
+    cwd = tmp_path / "out"
+    p = run_sos(repo_root / "pipeline/RNA_calling.ipynb", "rnaseqc_call",
+                {"sample_list": _two_sample_sheet(fix, tmp_path),
+                 "bam_list": tmp_path / "protocol_example.bam_file_list",
+                 "gtf": fix / "reference" / "chr22.window.collapsed.gtf",
+                 "reference_fasta": fix / "reference" / "chr22.window.fa.gz",
+                 "cwd": cwd, "numThreads": 1, "modular_script_dir": repo_root / "code/script"},
+                cwd=repo_root, timeout=600)
+    assert p.returncode == 0, p.stdout + p.stderr
+    for f in ("gene_tpm.gct.gz", "gene_readsCount.gct.gz", "exon_readsCount.gct.gz"):
+        with gzip.open(cwd / f"protocol_example.rnaseqc.{f}", "rt") as fh:
+            assert fh.readline().rstrip("\n").split("\t") == ["gene_ID", "SAMPLE_001", "SAMPLE_002"]
+    metrics = _rows((cwd / "protocol_example.rnaseqc.metrics.tsv").read_text())
+    assert [r[0] for r in metrics[1:]] == ["SAMPLE_001.cord.bam", "SAMPLE_002.cord.bam"]
diff --git a/tests/notebooks/pecotmr_integration/test_twas_ctwas.py b/tests/notebooks/pecotmr_integration/test_twas_ctwas.py
index 2d105f547..a68bcaeef 100644
--- a/tests/notebooks/pecotmr_integration/test_twas_ctwas.py
+++ b/tests/notebooks/pecotmr_integration/test_twas_ctwas.py
@@ -1,45 +1,96 @@
-"""Notebook tier: twas_ctwas.ipynb cTWAS chain on the blessed strong-signal fixture.
+"""Notebook tier: twas_ctwas.ipynb [twas] step on the toy TWAS weights, and the
+cTWAS chain on the blessed strong-signal fixture.
 
-Runs the full assemble -> estimate-priors -> fine-map chain over the 20-block
-chr22 LD grid with the blessed ctwas weights, reproducing the documented
+[twas] scores one region's genes and must reproduce the twas.R expected output.
+The cTWAS test runs the full assemble -> estimate-priors -> fine-map chain over
+the 20-block chr22 LD grid with the blessed ctwas weights, reproducing the documented
 gene-Z 3.168 / susie_pip 1.0. Exercises: the shared-var get_analysis_regions
 provider, the block-grid-vs-sketch collapse in gwas_sumstats_construct.R, the
 region provenance carried on the committed S4 TwasWeights fixture, and
 assembleCtwasInputs / estCtwasParam / finemapCtwasRegions via the ctwas_* wrappers.
+The assemble step alone also runs with the weights given as a list file and as a
+quoted wildcard, and rejects the two weight options given together.
 """
 from __future__ import annotations
 
 import pytest
 
+from helpers.expected import assert_matches_expected
 
-def test_ctwas(run_sos, read_rds, repo_root, tmp_path):
+
+def test_twas(run_sos, repo_root, tmp_path):
     fx = repo_root / "tests" / "fixtures"
-    cwd = tmp_path / "ctwas"
-    # The blessed cTWAS weights as committed S4 TwasWeights (region-provenance
-    # stamped -- chrom/p0/p1 -- which assembleCtwasInputs needs for LD-block
-    # placement; the legacy_ctwas_weights_to_s4 conversion step is retired).
-    s4 = fx / "twas/protocol_example.ctwas_weights.s4.chr22.rds"
-
-    # Full cTWAS chain: assemble + estimate priors + fine-map. The 20-block
-    # ld_meta doubles as the block grid (ctwas_manifest.R) and the LD sketch
-    # (gwas_sumstats_construct.R collapses it to the one chr22 panel).
+    # The committed xQTL meta tables point at legacy weights; [twas] reads the S4
+    # TwasWeights that twas.R's own test scores.
+    tw = fx / "twas/protocol_example.twas.reshaped_toy.chr22_ENSG00000130538.twas_weights.s4.rds"
+    meta = tmp_path / "xqtl_meta.tsv"
+    meta.write_text("#chr\tstart\tend\tregion_id\tTSS\toriginal_data\tcontexts\n"
+                    f"chr22\t10000000\t19000000\tENSG00000130538\t15528191\t{tw}\tbulk_rnaseq\n")
+    cwd = tmp_path / "out"
     p = run_sos(
-        repo_root / "pipeline/twas_ctwas.ipynb", "ctwas",
+        repo_root / "pipeline/twas_ctwas.ipynb", "twas",
         {
-            "run_param_est": True,
-            "run_finemapping": True,
             "cwd": cwd,
             "name": "protocol_example",
             "gwas_meta_data": fx / "twas/protocol_example.twas.gwas_meta.tsv",
-            "xqtl_meta_data": fx / "twas/protocol_example.twas.xqtl_meta.TAB.tsv",
-            "ld_meta_data": fx / "ld_reference/ld_meta_file.ctwas.tsv",
+            "xqtl_meta_data": meta,
+            "ld_meta_data": fx / "ld_reference/ld_meta_file.tsv",
             "region-name": "chr22_10000000_19000000",
-            "twas_weights": s4,
-            # The committed weights carry susie only; pecotmr resolves the
-            # method implicitly just for an 'ensemble' entry.
-            "ctwas_method": "susie",
+            # twas.R's default, under which its expected output was made (the
+            # notebook defaults to 0.05)
+            "rsq_pval_cutoff": "inf",
             "modular_script_dir": repo_root / "code/script",
         },
+        cwd=repo_root, timeout=600)
+    assert p.returncode == 0, p.stdout + p.stderr
+
+    # The genes reach twas.R as a per-gene table (a chromosome-wide region's
+    # lists are too long for the command line).
+    region = cwd / "twas" / "protocol_example.chr22_10000000_19000000"
+    table = (region.parent / f"{region.name}.twas_genes.tsv").read_text().splitlines()
+    assert table[0] == "twas_weights\tfine_mapping_result\toutput"
+    assert [row.split("\t")[0] for row in table[1:]] == [str(tw)]
+    assert_matches_expected(region.parent / f"{region.name}.ENSG00000130538.twas.rds",
+                            fx / "twas/expected/twas.chr22.rds",
+                            mode="tolerant", rtol=1e-6, atol=1e-8)
+
+
+# The blessed cTWAS weights as committed S4 TwasWeights (region-provenance
+# stamped -- chrom/p0/p1 -- which assembleCtwasInputs needs for LD-block
+# placement; the legacy_ctwas_weights_to_s4 conversion step is retired).
+CTWAS_WEIGHTS = "twas/protocol_example.ctwas_weights.s4.chr22.rds"
+
+
+def _ctwas_params(repo_root, cwd, **weights):
+    """The ctwas workflow's options, with the weights given as --twas-weights
+    or --twas-weights-list."""
+    fx = repo_root / "tests" / "fixtures"
+    # The 20-block ld_meta doubles as the block grid (ctwas_manifest.R) and the
+    # LD sketch (gwas_sumstats_construct.R collapses it to the one chr22 panel).
+    return {
+        "cwd": cwd,
+        "name": "protocol_example",
+        "gwas_meta_data": fx / "twas/protocol_example.twas.gwas_meta.tsv",
+        "xqtl_meta_data": fx / "twas/protocol_example.twas.xqtl_meta.TAB.tsv",
+        "ld_meta_data": fx / "ld_reference/ld_meta_file.ctwas.tsv",
+        "region-name": "chr22_10000000_19000000",
+        **weights,
+        # The committed weights carry susie only; pecotmr resolves the
+        # method implicitly just for an 'ensemble' entry.
+        "ctwas_method": "susie",
+        "modular_script_dir": repo_root / "code/script",
+    }
+
+
+def test_ctwas(run_sos, read_rds, repo_root, tmp_path):
+    cwd = tmp_path / "ctwas"
+    s4 = repo_root / "tests" / "fixtures" / CTWAS_WEIGHTS
+
+    # Full cTWAS chain: assemble + estimate priors + fine-map.
+    p = run_sos(
+        repo_root / "pipeline/twas_ctwas.ipynb", "ctwas",
+        {"run_param_est": True, "run_finemapping": True,
+         **_ctwas_params(repo_root, cwd, twas_weights=s4)},
         cwd=repo_root, timeout=1200)
     assert p.returncode == 0, p.stdout + p.stderr
 
@@ -53,3 +104,42 @@ def test_ctwas(run_sos, read_rds, repo_root, tmp_path):
     assert info["class"] == "CtwasResult"
     assert info["geneMaxPip"] == pytest.approx(1.0, abs=1e-3), info
     assert info["geneTopZ"] == pytest.approx(3.168, abs=0.01), info
+
+
+@pytest.mark.parametrize("source", ["twas_weights_list", "quoted_wildcard"])
+def test_ctwas_assemble_weights(run_sos, repo_root, tmp_path, source):
+    # A list file is how a genome-wide set of weights gets in (too long for the
+    # command line); a quoted wildcard is expanded by SoS. The folder name has a
+    # space, so the step must quote each path.
+    folder = tmp_path / "prebuilt weights"
+    folder.mkdir()
+    s4 = folder / "protocol_example.ctwas_weights.s4.chr22.rds"
+    s4.symlink_to(repo_root / "tests" / "fixtures" / CTWAS_WEIGHTS)
+    if source == "twas_weights_list":
+        listing = tmp_path / "weights.txt"
+        listing.write_text(f"{s4}\n")
+        weights = {"twas_weights_list": listing}
+    else:
+        weights = {"twas_weights": f"{folder}/*.rds"}
+    cwd = tmp_path / "ctwas"
+    p = run_sos(repo_root / "pipeline/twas_ctwas.ipynb", "ctwas",
+                _ctwas_params(repo_root, cwd, **weights),
+                cwd=repo_root, timeout=1200)
+    assert p.returncode == 0, p.stdout + p.stderr
+
+    # The step hands ctwas_assemble.R the same paths, one per line.
+    out = cwd / "ctwas"
+    assert (out / "protocol_example.ctwas_inputs.weights.txt").read_text().splitlines() == [str(s4)]
+    assert (out / "protocol_example.ctwas_inputs.rds").exists()
+
+
+def test_ctwas_weights_options_exclusive(run_sos, repo_root, tmp_path):
+    s4 = repo_root / "tests" / "fixtures" / CTWAS_WEIGHTS
+    listing = tmp_path / "weights.txt"
+    listing.write_text(f"{s4}\n")
+    p = run_sos(repo_root / "pipeline/twas_ctwas.ipynb", "ctwas",
+                _ctwas_params(repo_root, tmp_path / "ctwas",
+                              twas_weights=s4, twas_weights_list=listing),
+                cwd=repo_root, timeout=600)
+    assert p.returncode != 0
+    assert "Pass --twas-weights or --twas-weights-list, not both." in p.stdout + p.stderr
diff --git a/tests/scripts/pecotmr_integration/test_ctwas_assemble.py b/tests/scripts/pecotmr_integration/test_ctwas_assemble.py
index f90fcad72..c482b7db6 100644
--- a/tests/scripts/pecotmr_integration/test_ctwas_assemble.py
+++ b/tests/scripts/pecotmr_integration/test_ctwas_assemble.py
@@ -6,6 +6,7 @@
 from helpers.expected import assert_matches_expected
 
 EXPECTED = "tests/fixtures/ctwas/expected/ctwas_assemble.rds"
+WEIGHTS = "tests/fixtures/twas/protocol_example.ctwas_weights.s4.chr22.rds"
 
 
 def test_ctwas_assemble(ctwas_chain, read_rds, repo_root):
@@ -14,3 +15,20 @@ def test_ctwas_assemble(ctwas_chain, read_rds, repo_root):
     # embeds absolute LD/snp-map reference paths -> normalize_paths (basename compare).
     assert_matches_expected(ctwas_chain["inputs"], repo_root / EXPECTED,
                             mode="tolerant", rtol=1e-6, atol=1e-8, normalize_paths=True)
+
+
+def test_ctwas_assemble_list_file(ctwas_chain, run_r, repo_root, tmp_path):
+    """--twas-weights-list, the notebook's form (a genome-scale list is too long for
+    the command line), assembles the same inputs as --twas-weights; the arguments
+    otherwise match the ctwas_chain call."""
+    weights = tmp_path / "weights.txt"
+    weights.write_text(f"{repo_root / WEIGHTS}\n")
+    out = tmp_path / "ctwas_inputs.rds"
+    p = run_r(repo_root / "code/script/pecotmr_integration/ctwas_assemble.R",
+              ["--manifest", ctwas_chain["manifest"], "--twas-weights-list", weights,
+               "--method", "susie", "--twas-weight-cutoff", "0",
+               "--cs-min-cor", "0", "--min-pip-cutoff", "0", "--max-num-variants", "Inf",
+               "--output", out], timeout=300)
+    assert p.returncode == 0, p.stdout + p.stderr
+    assert_matches_expected(out, repo_root / EXPECTED,
+                            mode="tolerant", rtol=1e-6, atol=1e-8, normalize_paths=True)
diff --git a/tests/scripts/pecotmr_integration/test_fine_mapping_bed_merge.py b/tests/scripts/pecotmr_integration/test_fine_mapping_bed_merge.py
index 470675741..bd9a22cca 100644
--- a/tests/scripts/pecotmr_integration/test_fine_mapping_bed_merge.py
+++ b/tests/scripts/pecotmr_integration/test_fine_mapping_bed_merge.py
@@ -6,6 +6,8 @@
 import gzip
 import subprocess
 
+import pytest
+
 # The published vocabulary (pecotmr::exportToBed): the VCF's own field
 # names, shared by both schemas, not the retired variant_ID / event_ID /
 # region_ID / cs_coverage_0.95 spellings.
@@ -23,7 +25,8 @@ def _write_with(path, header, rows):
             handle.write("\t".join(str(v) for v in row) + "\n")
 
 
-def test_merge_sorts_numerically_and_indexes(run_r, repo_root, tmp_path):
+@pytest.mark.parametrize("via_list", [False, True], ids=["inline", "list_file"])
+def test_merge_sorts_numerically_and_indexes(via_list, run_r, repo_root, tmp_path):
     # Two regions, deliberately out of order and spanning chr2 vs chr10 so a
     # plain character sort would put 10 before 2.
     a = tmp_path / "r1.exported.bed.gz"
@@ -32,9 +35,16 @@ def test_merge_sorts_numerically_and_indexes(run_r, repo_root, tmp_path):
     _write(b, [("2", 100, 101, "2:101:A:T", "c1", "chr2_1_1000", 0.5),
                ("10", 50, 51, "10:51:A:T", "c1", "chr10_1_1000", 0.7)])
     out = tmp_path / "merged.bed.gz"
+    if via_list:
+        # the notebook's form: a genome-wide list is too long for the command line
+        listing = tmp_path / "inputs.txt"
+        listing.write_text(f"{a}\n{b}\n")
+        inputs = ["--input-list", listing]
+    else:
+        inputs = ["--input", a, b]
 
     p = run_r(repo_root / "code/script/pecotmr_integration/fine_mapping_bed_merge.R",
-              ["--input", a, b, "--output", out], timeout=300)
+              [*inputs, "--output", out], timeout=300)
     assert p.returncode == 0, p.stdout + p.stderr
 
     lines = gzip.open(out, "rt").read().rstrip("\n").split("\n")
@@ -127,3 +137,17 @@ def test_merge_does_not_index_a_position_table(run_r, repo_root, tmp_path):
     # alleles survive untouched in the allele column
     assert [r[2] for r in rows] == ["G", "A", "C"]
     assert not (tmp_path / "merged.positions.tsv.gz.tbi").exists()
+
+
+def test_merge_rejects_tables_with_different_headers(run_r, repo_root, tmp_path):
+    """The tables are read one at a time, so the merge checks that they share one
+    header: binding an fSuSiE table to a top-loci one (same width, different
+    columns) would otherwise mislabel or NA-pad the merged table."""
+    a = tmp_path / "r1.exported.bed.gz"
+    b = tmp_path / "r2.top_loci.bed.gz"
+    _write(a, [("2", 300, 301, "2:301:A:T", "c1", "chr2_1_1000", 0.9)])
+    _write_with(b, TOPLOCI_HEADER, [("chr2", 99, 100, "chr2:100:G:A", "G", "A", 0.5)])
+    p = run_r(repo_root / "code/script/pecotmr_integration/fine_mapping_bed_merge.R",
+              ["--input", a, b, "--output", tmp_path / "merged.bed.gz"], timeout=300)
+    assert p.returncode != 0
+    assert "must share the header" in p.stderr, p.stderr
diff --git a/tests/scripts/pecotmr_integration/test_fine_mapping_upset.py b/tests/scripts/pecotmr_integration/test_fine_mapping_upset.py
index 4cb4c1627..85bfaed2c 100644
--- a/tests/scripts/pecotmr_integration/test_fine_mapping_upset.py
+++ b/tests/scripts/pecotmr_integration/test_fine_mapping_upset.py
@@ -12,10 +12,18 @@
 import pytest
 
 
-def test_fine_mapping_upset(fmr, run_r, repo_root, tmp_path):
+@pytest.mark.parametrize("via_list", [False, True], ids=["inline", "list_file"])
+def test_fine_mapping_upset(via_list, fmr, run_r, repo_root, tmp_path):
     out = tmp_path / "upset.png"
+    if via_list:
+        # the notebook's form: a genome-wide list is too long for the command line
+        listing = tmp_path / "inputs.txt"
+        listing.write_text(f"{fmr}\n")
+        inputs = ["--input-list", listing]
+    else:
+        inputs = ["--input", fmr]
     p = run_r(repo_root / "code/script/pecotmr_integration/fine_mapping_upset.R",
-              ["--input", fmr, "--output", out], timeout=120)
+              [*inputs, "--output", out], timeout=120)
     assert p.returncode == 0, p.stdout + p.stderr
     assert out.exists() and out.stat().st_size > 0
     assert "No credible sets found across any input" in p.stderr, (
diff --git a/tests/scripts/pecotmr_integration/test_mash_feature_score.py b/tests/scripts/pecotmr_integration/test_mash_feature_score.py
index e7d6e3af5..093beaae1 100644
--- a/tests/scripts/pecotmr_integration/test_mash_feature_score.py
+++ b/tests/scripts/pecotmr_integration/test_mash_feature_score.py
@@ -54,7 +54,8 @@ def test_mash_feature_score(method, score_type, run_r, read_rds, repo_root, tmp_
     # write.csv(row.names=FALSE) drops those path-derived names.
 
 
-def test_mash_feature_score_merge(run_r, repo_root, tmp_path):
+@pytest.mark.parametrize("via_list", [False, True], ids=["inline", "list_file"])
+def test_mash_feature_score_merge(via_list, run_r, repo_root, tmp_path):
     contrast = _contrast(run_r, repo_root, tmp_path)
     fs = tmp_path / "fs.rds"
     p = run_r(repo_root / "code/script/pecotmr_integration/mash_feature_score.R",
@@ -63,8 +64,15 @@ def test_mash_feature_score_merge(run_r, repo_root, tmp_path):
     assert p.returncode == 0, p.stdout + p.stderr
 
     out = tmp_path / "merged.csv"
+    if via_list:
+        # the notebook's form: a genome-wide list is too long for the command line
+        listing = tmp_path / "scores.txt"
+        listing.write_text(f"{fs}\n")
+        scores = ["--scores-list", listing]
+    else:
+        scores = ["--scores", fs]
     p = run_r(repo_root / "code/script/pecotmr_integration/mash_feature_score_merge.R",
-              ["--scores", fs, "--output", out], timeout=120)
+              [*scores, "--output", out], timeout=120)
     assert p.returncode == 0, p.stdout + p.stderr
     assert out.exists()
     assert "scoreType" in out.read_text().splitlines()[0]
diff --git a/tests/scripts/pecotmr_integration/test_mash_posterior_contrast.py b/tests/scripts/pecotmr_integration/test_mash_posterior_contrast.py
index a535d820d..e525d8c97 100644
--- a/tests/scripts/pecotmr_integration/test_mash_posterior_contrast.py
+++ b/tests/scripts/pecotmr_integration/test_mash_posterior_contrast.py
@@ -55,3 +55,24 @@ def test_mash_posterior_contrast_summary_and_plot(run_r, repo_root, tmp_path):
               ["--data", summ, "--output", png], timeout=300)
     assert p.returncode == 0, p.stdout + p.stderr
     assert png.exists() and png.stat().st_size > 1000
+
+
+def test_mash_posterior_contrast_summary_list_file(run_r, repo_root, tmp_path):
+    """--contrast-list (the notebook's form for a genome-wide list, too long for the
+    command line) tallies the same regions as --contrast; blank lines are skipped."""
+    contrast = _contrast(run_r, repo_root, tmp_path)
+    script = repo_root / "code/script/pecotmr_integration/mash_posterior_contrast_summary.R"
+    inline = tmp_path / "inline.csv"
+    p = run_r(script, ["--contrast", contrast, contrast, "--cells", CELLS,
+                       "--p-cutoff", "0.5", "--output", inline], timeout=300)
+    assert p.returncode == 0, p.stdout + p.stderr
+    listing = tmp_path / "contrast.txt"
+    listing.write_text(f"{contrast}\n\n{contrast}\n")
+    listed = tmp_path / "listed.csv"
+    p = run_r(script, ["--contrast-list", listing, "--cells", CELLS,
+                       "--p-cutoff", "0.5", "--output", listed], timeout=300)
+    assert p.returncode == 0, p.stdout + p.stderr
+    assert listed.read_text() == inline.read_text()
+    n_all = next(line for line in listed.read_text().splitlines()
+                 if line.startswith('"n_all_feature"'))
+    assert set(n_all.split(",")[1:]) == {"2"}, "both listed regions are counted"
diff --git a/tests/scripts/pecotmr_integration/test_mash_preprocessing.py b/tests/scripts/pecotmr_integration/test_mash_preprocessing.py
index 545295e13..85cfd6d51 100644
--- a/tests/scripts/pecotmr_integration/test_mash_preprocessing.py
+++ b/tests/scripts/pecotmr_integration/test_mash_preprocessing.py
@@ -13,6 +13,8 @@
 """
 from __future__ import annotations
 
+import pytest
+
 from helpers.expected import assert_matches_expected
 
 MASH_KEYS = {
@@ -49,11 +51,23 @@ def test_mash_preprocessing_qtlsumstats_path(run_r, read_rds, repo_root, tmp_pat
                             rtol=1e-6, atol=1e-8)
 
 
-def test_mash_preprocessing_finemapping_path(run_r, read_rds, repo_root, tmp_path):
+@pytest.mark.parametrize("via_list", [False, True], ids=["inline", "list_file"])
+def test_mash_preprocessing_finemapping_path(via_list, run_r, read_rds, repo_root, tmp_path):
     fx = repo_root / "tests/fixtures/mash"
+    fmr = fx / "protocol_example.QtlFineMappingResult.rds"
     out = tmp_path / "mash_input_fmr.rds"
+    if via_list:
+        # the notebook's form for a genome-wide list. The label is the default
+        # (the file basename), so the result matches the inline expected output.
+        objects = tmp_path / "objects.txt"
+        objects.write_text(f"{fmr}\n")
+        labels = tmp_path / "region_ids.txt"
+        labels.write_text("protocol_example.QtlFineMappingResult\n")
+        inputs = ["--objects-list", objects, "--region-ids-list", labels]
+    else:
+        inputs = ["--objects", fmr]
     p = run_r(repo_root / "code/script/pecotmr_integration/mash_preprocessing.R",
-              ["--objects", fx / "protocol_example.QtlFineMappingResult.rds",
+              [*inputs,
                "--n-random", "15", "--n-null", "15", "--coverage", "0.95",
                "--sig-p-cutoff", "0.1", "--seed", "1", "--output", out], timeout=300)
     assert p.returncode == 0, p.stdout + p.stderr
diff --git a/tests/scripts/pecotmr_integration/test_qtl_enrichment.py b/tests/scripts/pecotmr_integration/test_qtl_enrichment.py
index 84fe3a389..2fe8f2566 100644
--- a/tests/scripts/pecotmr_integration/test_qtl_enrichment.py
+++ b/tests/scripts/pecotmr_integration/test_qtl_enrichment.py
@@ -2,16 +2,28 @@
 QTL + GWAS S4 FineMappingResult fixtures (legacy conversion retired)."""
 from __future__ import annotations
 
+import pytest
+
 FX = "tests/fixtures/susie_enloc"
 QTL = FX + "/protocol_example.enloc.MiGA_eQTL.ENSG00000142798.fine_mapping.rds"
 GWAS = FX + "/protocol_example.enloc.gwas.chr1_20110062_22020160.fine_mapping.rds"
 
 
-def test_qtl_enrichment(run_r, read_rds, repo_root, tmp_path):
+@pytest.mark.parametrize("via_list", [False, True], ids=["inline", "list_file"])
+def test_qtl_enrichment(via_list, run_r, read_rds, repo_root, tmp_path):
     out = tmp_path / "enrich.rds"
+    if via_list:
+        # the notebook's form: a condition's genome-wide list is too long for the
+        # command line
+        qtl_list = tmp_path / "qtl.txt"
+        qtl_list.write_text(f"{repo_root / QTL}\n")
+        gwas_list = tmp_path / "gwas.txt"
+        gwas_list.write_text(f"{repo_root / GWAS}\n")
+        inputs = ["--qtl-fine-mapping-list", qtl_list, "--gwas-fine-mapping-list", gwas_list]
+    else:
+        inputs = ["--qtl-fine-mapping", repo_root / QTL, "--gwas-fine-mapping", repo_root / GWAS]
     p = run_r(repo_root / "code/script/pecotmr_integration/qtl_enrichment.R",
-              ["--qtl-fine-mapping", repo_root / QTL,
-               "--gwas-fine-mapping", repo_root / GWAS,
+              [*inputs,
                # --seed exercises the forwarding into qtlEnrichmentPipeline(seed=);
                # the C++ imputation sampler is unreachable from R's set.seed, so a
                # rejected/renamed argument would only ever surface here. Not a value
diff --git a/tests/scripts/pecotmr_integration/test_region_manifest.py b/tests/scripts/pecotmr_integration/test_region_manifest.py
index 279e0862e..78f2ab2c8 100644
--- a/tests/scripts/pecotmr_integration/test_region_manifest.py
+++ b/tests/scripts/pecotmr_integration/test_region_manifest.py
@@ -6,11 +6,12 @@
 from helpers.expected import assert_matches_expected
 
 EXPECTED = "tests/fixtures/region_manifest/expected/region_manifest.tsv"
+SCRIPT = "code/script/pecotmr_integration/region_manifest.R"
 
 
 def test_region_manifest(run_r, repo_root, qtl_mini, tmp_path):
     out = tmp_path / "regions.tsv"
-    p = run_r(repo_root / "code/script/pecotmr_integration/region_manifest.R",
+    p = run_r(repo_root / SCRIPT,
               ["--pheno-manifest", qtl_mini / "protocol_example.pheno_manifest_context.tsv",
                "--customized-association-windows", qtl_mini / "association_windows.bed",
                "--output", out], timeout=120)
@@ -21,3 +22,67 @@ def test_region_manifest(run_r, repo_root, qtl_mini, tmp_path):
     assert {"region_id", "chr", "start", "end"} <= set(cols)
     assert any("ENSG00000283047" in ln for ln in lines[1:])
     assert_matches_expected(out, repo_root / EXPECTED, mode="tolerant", rtol=1e-6, atol=1e-8)
+
+
+def _manifest_rows(run_r, repo_root, tmp_path, pheno, *args):
+    out = tmp_path / "manifest.tsv"
+    p = run_r(repo_root / SCRIPT,
+              ["--pheno-manifest", pheno, *args, "--output", out], timeout=120)
+    assert p.returncode == 0, p.stdout + p.stderr
+    rows = [line.split("\t") for line in out.read_text().splitlines()]
+    return [dict(zip(rows[0], row)) for row in rows[1:]]
+
+
+def test_region_manifest_selection(run_r, repo_root, tmp_path):
+    """Gene mode ([susie_twas] / [mnm]) runs on region_name, so it must keep an
+    ID as given (a versioned ENSG ID) while region_id is made filename-safe.
+    --region-name-list reads like --region-list and, like --region-name, wins
+    over it."""
+    pheno = tmp_path / "pheno.tsv"
+    pheno.write_text("#chr\tstart\tend\tID\n"
+                     "chr22\t100\t200\tENSG00000000001.5\n"
+                     "chr22\t300\t400\tENSG00000000002\n")
+    names = tmp_path / "names.txt"
+    names.write_text("# genes to run\n\nchr22 0 0 ENSG00000000002\n"
+                     "  ENSG00000000001.5  \nENSG00000000002\n")
+    regions = tmp_path / "regions.txt"
+    regions.write_text("ENSG00000000001.5\n")
+
+    rows = _manifest_rows(run_r, repo_root, tmp_path, pheno,
+                          "--region-name-list", names, "--region-list", regions)
+    assert [r["region_name"] for r in rows] == ["ENSG00000000002", "ENSG00000000001.5"]
+    assert [r["region_id"] for r in rows] == ["ENSG00000000002", "ENSG00000000001_5"]
+    rows = _manifest_rows(run_r, repo_root, tmp_path, pheno, "--region-list", regions)
+    assert [r["region_name"] for r in rows] == ["ENSG00000000001.5"]
+
+
+def test_region_manifest_gene_mode_without_coordinates(run_r, repo_root, tmp_path):
+    """qtl_dataset_construct does not need #chr/start/end, and neither does gene
+    mode: it lists the genes with empty coordinates. Region mode still needs them
+    for each gene's window, and gene mode still rejects an ID that is not a gene."""
+    pheno = tmp_path / "pheno.tsv"
+    pheno.write_text("ID\tpath\tcond\n"
+                     "ENSG00000000001.5\tx.bed.gz\tcontext1\n"
+                     "ENSG00000000002\tx.bed.gz\tcontext1\n")
+    rows = _manifest_rows(run_r, repo_root, tmp_path, pheno, "--gene-mode")
+    assert [r["region_name"] for r in rows] == ["ENSG00000000001.5", "ENSG00000000002"]
+    assert all(r["chr"] == r["start"] == r["end"] == r["ld_block"] == "" for r in rows)
+    for args, error in [((), "needs ID, #chr/chrom, start, end columns"),
+                        (("--gene-mode", "--region-name", "ENSG00000000009"),
+                         "Gene 'ENSG00000000009' not found in --pheno-manifest")]:
+        p = run_r(repo_root / SCRIPT, ["--pheno-manifest", pheno, *args,
+                                       "--output", tmp_path / "rejected.tsv"], timeout=120)
+        assert p.returncode != 0
+        assert error in p.stdout + p.stderr
+
+
+@pytest.mark.parametrize("flag", ["--region-list", "--customized-association-windows"])
+def test_region_manifest_missing_file(run_r, repo_root, qtl_mini, tmp_path, flag):
+    """A list or windows file that is given must exist: skipping it would
+    quietly analyze every region instead."""
+    p = run_r(repo_root / SCRIPT,
+              ["--pheno-manifest", qtl_mini / "protocol_example.pheno_manifest_context.tsv",
+               flag, tmp_path / "missing.txt", "--output", tmp_path / "manifest.tsv"],
+              timeout=120)
+    assert p.returncode != 0
+    assert "not found" in p.stdout + p.stderr
diff --git a/tests/scripts/pecotmr_integration/test_twas.py b/tests/scripts/pecotmr_integration/test_twas.py
index b8c7059ef..54d74c18a 100644
--- a/tests/scripts/pecotmr_integration/test_twas.py
+++ b/tests/scripts/pecotmr_integration/test_twas.py
@@ -2,6 +2,8 @@
 GwasSumStats built here (the legacy twas-weights conversion has been retired)."""
 from __future__ import annotations
 
+import pytest
+
 from helpers.expected import assert_matches_expected
 
 FX = "tests/fixtures"
@@ -9,7 +11,8 @@
 EXPECTED = FX + "/twas/expected/twas.chr22.rds"
 
 
-def test_twas(run_r, read_rds, repo_root, tmp_path):
+@pytest.mark.parametrize("via_table", [False, True], ids=["inline", "gene_table"])
+def test_twas(via_table, run_r, read_rds, repo_root, tmp_path):
     r = repo_root
     gss = tmp_path / "gss.rds"
     assert run_r(r / "code/script/pecotmr_integration/gwas_sumstats_construct.R",
@@ -18,8 +21,16 @@ def test_twas(run_r, read_rds, repo_root, tmp_path):
          "--ld-block", "chr22:10000000-19000000",
          "--ld-meta", r / f"{FX}/ld_reference/ld_meta_file.tsv", "--output", gss], timeout=300).returncode == 0
     out = tmp_path / "twas.rds"
+    if via_table:
+        # the notebook's form: a chromosome-wide region's per-gene lists are too
+        # long for the command line. An empty fine_mapping_result means no MR.
+        table = tmp_path / "genes.tsv"
+        table.write_text(f"twas_weights\tfine_mapping_result\toutput\n{r / TW}\t\t{out}\n")
+        genes = ["--gene-table", table]
+    else:
+        genes = ["--twas-weights", r / TW, "--output", out]
     p = run_r(r / "code/script/pecotmr_integration/twas.R",
-              ["--twas-weights", r / TW, "--gwas-sumstats", gss, "--output", out], timeout=300)
+              [*genes, "--gwas-sumstats", gss], timeout=300)
     assert p.returncode == 0, p.stdout + p.stderr
     assert read_rds(out)["class"] == "GRanges"
     assert_matches_expected(out, repo_root / EXPECTED, mode="tolerant", rtol=1e-6, atol=1e-8)