Feat/pipeline selector landing page

code/SoS/association_scan/qtl_association_postprocessing.ipynb

@@ -8,7 +8,7 @@
    "source": [
     "# Hierarchical Multiple Testing\n",
     "\n",
-    "This is a runnable, toy-data example of the `qtl_association_testing` post-processing step. It performs hierarchical multiple-testing correction on the per-gene cis-QTL association results produced by TensorQTL, then archives the bulky intermediate files."
+    "Applies hierarchical multiple-testing correction to the per-gene cis-QTL results from TensorQTL and assembles the survivors into one consolidated object."
    ]
   },
   {
@@ -17,11 +17,25 @@
     "kernel": "SoS"
    },
    "source": [
-    "## Description\n",
+    "## Overview\n",
     "\n",
-    "Given the per-gene association output from a cis-QTL scan, this pipeline implements a three-step hierarchical procedure: (1) local adjustment of the p-values of all cis-variants within each gene, (2) global adjustment of the minimum adjusted p-value across genes, and (3) identification of significant xQTL whose locally adjusted p-value falls below the threshold. It also reorganises the intermediate TensorQTL files into an archive folder for book-keeping or deletion.\n",
+    "A cis scan reports a p-value for every variant against every molecular phenotype, which is\n",
+    "not yet a result: the variants within a gene's window are correlated and the genes are\n",
+    "many, so raw p-values overstate significance twice over.\n",
     "\n",
-    "The same `default` workflow handles three flavours of upstream results: the **genetic (cis) effect** (shown below), the **interaction effect** (run on `*.cis_qtl_top_assoc.txt.gz` with `--maf-cutoff 0 --cis-window 0` and interaction p/q-value patterns), and **quantile QTL**. This example demonstrates the genetic-effect flavour on toy data. (cf. Huang et al. 2018)"
+    "The correction is hierarchical, in three steps: local adjustment of the p-values of all cis\n",
+    "variants within each gene, global adjustment of the minimum adjusted p-value across genes,\n",
+    "then selection of the xQTLs whose locally adjusted p-value falls below the threshold\n",
+    "(the eigenMT-BH procedure of [Huang et al. 2018, *NAR* 46(22):e133](https://doi.org/10.1093/nar/gky780)). The survivors are packaged as a `QtlSumStats` object with a\n",
+    "regional FDR table, and the intermediate TensorQTL files are reorganised into an archive\n",
+    "folder for book-keeping or deletion.\n",
+    "\n",
+    "The same `default` workflow handles three flavours of upstream result: the genetic (cis)\n",
+    "effect shown below, the interaction effect - run on `*.cis_qtl_top_assoc.txt.gz` with\n",
+    "`--maf-cutoff 0 --cis-window 0` and the interaction p/q-value patterns - and quantile QTL.\n",
+    "\n",
+    "**When to run it.** Immediately after `TensorQTL`, before fine-mapping. Nothing downstream\n",
+    "reads the raw scan directory directly."
    ]
   },
   {
@@ -30,7 +44,44 @@
     "kernel": "SoS"
    },
    "source": [
-    "**Timing:** ~1-2 min on typical compute infrastructure."
+    "## Input\n",
+    "\n",
+    "The workflow is pointed at a directory rather than at files: it globs `--cwd` for the three TensorQTL products below using `--regional-pattern`, `--qtl-pattern` and `--n-variants-suffix`. `output/tensorqtl_cis` is one directory in that shape.\n",
+    "\n",
+    "- `--cwd`: the TensorQTL working directory holding the scan output. Three kinds of file are picked up from it.\n",
+    "\n",
+    "  The per-gene regional result, matched by `--regional-pattern` (default `cis_qtl[.]regional[.]tsv[.]gz$`). Example `output/tensorqtl_cis/protocol_example.chr22_chr22.cis_qtl.regional.tsv.gz`, 21 columns and 201 rows:\n",
+    "\n",
+    "  ```\n",
+    "  chrom  pos       n_variants  beta_shape1  beta_shape2  true_df    p_true_df             variant_id\n",
+    "  22     11797697  6           0.9833031    5.6122484    30.858135  0.0676042590649028    chr22:11797697_C_T\n",
+    "  22     15838003  46          0.98497593   26.038029    30.771364  0.007232871450238482  chr22:15838003_CTTTC...\n",
+    "  ```\n",
+    "\n",
+    "  The full variant-level pairs file, matched by `--qtl-pattern` (default `cis_qtl[.]pairs[.]tsv[.]gz$`). Example `output/tensorqtl_cis/protocol_example.chr22_chr22.cis_qtl.pairs.tsv.gz`, 15 columns and 18319 rows:\n",
+    "\n",
+    "  ```\n",
+    "  chrom  pos       molecular_trait_id  variant_id          tss_distance  tes_distance  af           pvalue\n",
+    "  22     10685239  ENSG00000283047     chr22:10685239_C_T  -254148       -276098       0.059322033  0.572662489372593\n",
+    "  22     10761227  ENSG00000283047     chr22:10761227_G_A  -178160       -200110       0.15254237   0.752156227936135\n",
+    "  ```\n",
+    "\n",
+    "  The recounted variant tallies, matched by `--n-variants-suffix` (default `cis_n_variants_stats[.]tsv[.]gz$`). Example `output/tensorqtl_cis/protocol_example.chr22_chr22.maf_0.01_window_1000000_cis_n_variants_stats.tsv.gz`, 4 columns and 201 rows:\n",
+    "\n",
+    "  ```\n",
+    "  chrom  molecular_trait_object_id  n_variants  n_variants_filtered\n",
+    "  chr22  ENSG00000063515            44          44\n",
+    "  chr22  ENSG00000075234            111         111\n",
+    "  ```\n",
+    "\n",
+    "- `--output-dir`: where the consolidated object and the regional FDR table are written. Required.\n",
+    "- `--modular-script-dir`: the `code/script` directory holding the pecotmr_integration wrappers. Required.\n",
+    "- `--maf-cutoff` and `--cis-window`: `0.01` and `1000000` by default. When either is non-zero the variant counts are recomputed from the full pairs file before Bonferroni correction, which matters because a correction taken over too wide a window, or over many low-frequency variants, overstates the number of independent tests. Set both to `0` to skip the recount and use the counts as given.\n",
+    "- `--pvalue-cutoff`: `0.05` by default. Pairs above it are dropped before assembly to keep the object small; the per-gene minimum variant is always retained so the Bonferroni minimum stays exact.\n",
+    "- `--fdr-threshold`: `0.05` by default, the threshold for the global significance call.\n",
+    "- `--study`, `--context`, `--genome`: labels stamped onto the assembled `QtlSumStats` object, `study` / `context` / `hg38` by default.\n",
+    "- `--af-col`, `--pvalue-col`, `--molecular-id-col`: the column names to read, `af` / `pvalue` / `molecular_trait_object_id` by default.\n",
+    "- `--[no-]enable-archive`: off by default. When on, the bulky TensorQTL intermediates are moved out of the working directory for book-keeping or deletion."
    ]
   },
   {
@@ -39,13 +90,14 @@
     "kernel": "SoS"
    },
    "source": [
-    "## Input Files\n",
+    "## Output\n",
+    "\n",
+    "- `{output-dir}/{cwd basename}.qtl_association_postprocessing.rds` - the consolidated `QtlSumStats` object holding the surviving xQTLs together with the `--study`, `--context` and `--genome` labels. This is the only file named in the workflow's `output:` statement.\n",
+    "- `{output-dir}/{cwd basename}.qtl_association_postprocessing.cis_regional.fdr.tsv.gz` - the per-gene regional FDR table.\n",
+    "\n",
+    "Both are named from the basename of `--cwd`, so the example below writes `tensorqtl_cis.qtl_association_postprocessing.rds` and `tensorqtl_cis.qtl_association_postprocessing.cis_regional.fdr.tsv.gz` into `output/qtlpp`.\n",
     "\n",
-    "| Role | Toy file | Produced by |\n",
-    "|------|----------|-------------|\n",
-    "| Per-gene cis association | `output/tensorqtl_cis/*.cis_qtl.regional.tsv.gz` | TensorQTL cis scan |\n",
-    "| Full variant-level pairs | `output/tensorqtl_cis/*.cis_qtl.pairs.tsv.gz` | TensorQTL cis scan |\n",
-    "| Gene coordinates | `input/reference_data/look_up_gene_id.tsv` (`chr`, `start`, `end`, `gene_id`) | reference |"
+    "No example output ships for this step. The `*.cis_regional.*.fdr.gz` and `*_regional_significance.*` files sitting alongside the inputs in `output/tensorqtl_cis` were written by an earlier version of this interface and carry names the workflow no longer produces."
    ]
   },
   {
@@ -54,19 +106,9 @@
     "kernel": "SoS"
    },
    "source": [
-    "## Output Files\n",
+    "## Minimal Working Example\n",
     "\n",
-    "Written under `--output-dir` (significance results) and `--archive-dir` (moved intermediates):\n",
-    "\n",
-    "| File | Meaning |\n",
-    "|------|---------|\n",
-    "| `*_multiple_testing_consolidated.rds` | Consolidated multiple-testing results object |\n",
-    "| `*.cis_regional.summary.txt` | Per-gene regional summary after correction |\n",
-    "| `*significant_qtl*.tsv.gz` | Significant variant-level QTL (Bonferroni / BH / permutation) |\n",
-    "| `*significant_events*.tsv.gz` | Significant gene-level events |\n",
-    "| `archive/.../*.parquet` | Bulky TensorQTL intermediates moved out of the working dir |\n",
-    "\n",
-    "On the 6-sample toy data the correction runs cleanly and produces the consolidated object and summary; no events pass the genome-wide significance threshold, which is expected at this tiny scale."
+    "The workflow globs `--cwd` for its inputs rather than taking file arguments, so it is pointed straight at the TensorQTL scan directory."
    ]
   },
   {
@@ -75,9 +117,7 @@
     "kernel": "SoS"
    },
    "source": [
-    "## Steps\n",
-    "\n",
-    "Run the hierarchical correction on the toy cis results. The MAF / cis-window options are non-zero, so the program recomputes the variant counts from the full pairs file before applying the Bonferroni correction. On this toy data the distance columns in the TensorQTL output are named `tss_distance` and `tes_distance`, so we point `--tss-dist-col` / `--tes-dist-col` at them."
+    "**Timing**: TBD (on toy dataset)"
    ]
   },
   {
@@ -90,15 +130,10 @@
    "source": [
     "sos run pipeline/qtl_association_postprocessing.ipynb default \\\n",
     "    --cwd output/tensorqtl_cis \\\n",
-    "    --gene-coordinates input/reference_data/look_up_gene_id.tsv \\\n",
-    "    --sub-dir . \\\n",
-    "    --tss-dist-col tss_distance --tes-dist-col tes_distance \\\n",
-    "    --maf-cutoff 0.01 --cis-window 1000000 \\\n",
-    "    --regional-pattern \"*.cis_qtl.regional.tsv.gz$\" \\\n",
-    "    --output-dir output/hierarchical_multi_test/output \\\n",
-    "    --archive-dir output/hierarchical_multi_test/archive --enable-archive True \\\n",
-    "    --pecotmr-path ../pecotmr \\\n",
-    "    -s force"
+    "    --output-dir output/qtlpp \\\n",
+    "    --modular-script-dir code/script \\\n",
+    "    --maf-cutoff 0.01 --cis-window 1000000 --pvalue-cutoff 0.05 \\\n",
+    "    --study protocol_example --context bulk_rnaseq --genome hg38"
    ]
   },
   {
@@ -107,16 +142,14 @@
     "kernel": "SoS"
    },
    "source": [
-    "## Command interface\n",
-    "\n",
-    "Show all workflow options and their defaults:"
+    "## Command Interface"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {
-    "kernel": "Bash"
+    "kernel": "SoS"
    },
    "outputs": [],
    "source": [
@@ -129,9 +162,48 @@
     "kernel": "SoS"
    },
    "source": [
-    "## Workflow implementation\n",
+    "```\n",
+    "usage: sos run pipeline/qtl_association_postprocessing.ipynb\n",
+    "               [workflow_name | -t targets] [options] [workflow_options]\n",
+    "  workflow_name:        Single or combined workflows defined in this script\n",
+    "  targets:              One or more targets to generate\n",
+    "  options:              Single-hyphen sos parameters (see \"sos run -h\" for details)\n",
+    "  workflow_options:     Double-hyphen workflow-specific parameters\n",
     "\n",
-    "The cells below are the unmodified SoS workflow definition that the example above calls."
+    "Workflows:\n",
+    "  default\n",
+    "\n",
+    "Global Workflow Options:\n",
+    "  --cwd . (as path)\n",
+    "  --modular-script-dir VAL (as path, required)\n",
+    "                        code/script dir that holds the pecotmr_integration\n",
+    "                        wrappers\n",
+    "  --output-dir VAL (as path, required)\n",
+    "  --sub-dir . (as path)\n",
+    "  --study study\n",
+    "                        labels for the assembled QtlSumStats\n",
+    "  --context context\n",
+    "  --genome hg38\n",
+    "  --maf-cutoff 0.01 (as float)\n",
+    "  --cis-window 1000000 (as int)\n",
+    "  --af-col af\n",
+    "  --pvalue-col pvalue\n",
+    "  --molecular-id-col 'molecular_trait_object_id'\n",
+    "  --pvalue-cutoff 0.05 (as float)\n",
+    "                        drop pairs with p above this before assembling (keeps\n",
+    "                        the object small; the per-gene min variant is always\n",
+    "                        retained so the Bonferroni min is exact)\n",
+    "  --fdr-threshold 0.05 (as float)\n",
+    "  --regional-pattern 'cis_qtl[.]regional[.]tsv[.]gz$'\n",
+    "                        list.files() regex patterns to locate the inputs in the\n",
+    "                        work dir\n",
+    "  --qtl-pattern 'cis_qtl[.]pairs[.]tsv[.]gz$'\n",
+    "  --n-variants-suffix 'cis_n_variants_stats[.]tsv[.]gz$'\n",
+    "  --[no-]enable-archive (default to False)\n",
+    "\n",
+    "Sections\n",
+    "  default:\n",
+    "```"
    ]
   },
   {
@@ -140,9 +212,9 @@
     "kernel": "SoS"
    },
    "source": [
-    "## Anticipated Results\n",
+    "## Workflow implementation\n",
     "\n",
-    "The pipeline produces output files in the `output/` subdirectory named after the workflow step. Verify success by checking that output files exist and are non-empty. See the **Output** section above for the expected file names and formats."
+    "The cells below are the unmodified SoS workflow definition that the example above calls."
    ]
   },
   {