UC3 Setup Recipe — Differential ATAC-seq Accessibility

Issue: https://github.com/jmchilton/galaxy-brain/issues/14 Science: Genome-wide differential ATAC-seq accessibility, Erythroblast vs B-cell (Corces 2016, hg19). 45,620 significant peaks; GATA1/HBB Erythroblast-gained; PAX5/EBF1/MS4A1/CD19 B-cell-gained. Extracted workflow: 13 steps — counts collection + sample sheet → DESeq2 → NA-filter → { volcano ∥ significance-filter → rank → top-gained / top-lost }, with each PDF figure (PCA, volcano) converted to PNG in-graph (graphicsmagick → extract-page). 6 workflow outputs; re-runs to the identical numbers. Latest extraction: UC3_ATAC_extracted_figures.ga.

Server-agnostic. Tool versions confirmed on Galaxy 26.2.dev0. Build everything through /api/tools (or the MCP run_tool); the notebook and the extraction are plain API calls — no UI needed. No Galaxy code change required — DESeq2/volcano emit PDF, so each figure is rasterized to PNG by an in-graph tool (§1) and referenced with the stock history_dataset_as_image directive.

1. Figure rendering — PDF → image via in-graph tool (no code change)

The DESeq2 diagnostics and the volcano are PDF, which the notebook cannot embed as an image. Instead of a core PDF-rasterizing directive (the abandoned history_dataset_as_pdf / PyMuPDF path — see PDF_IMAGES.md), rasterize each figure with a tool, so the image becomes a real on-graph dataset referenced by the stock history_dataset_as_image directive. Provenance-clean, extraction-automatic, no server dependency.

Per PDF figure:

1. graphicsmagick_image_convert (bgruening) on the PDF, output_format=png. For PDF input it splits every page into a list collection splitted_pdf, elements temp_000, temp_001, … (one per page). Runs in the graphicsmagick biocontainer (ghostscript delegate reads the PDF).
2. __EXTRACT_DATASET__ by identifier to pull the wanted page (temp_000 = page 1) out of the collection as a standalone PNG HDA.
3. Rename that HDA to a meaningful figure name (Sample PCA plot, Volcano plot) — it becomes the workflow output label and the output="…" reference in the extracted report.

Volcano PDF = 1 page → temp_000. DESeq2 diagnostics = 5 pages → page 1 (temp_000) is the PCA. (Cost: graphicsmagick converts all pages, so the 5-page DESeq2 PDF yields 5 PNGs even though only the PCA is used.)

2. Tools

3. Input data

Source: Corces 2016 hematopoiesis ATAC count matrix — GEO GSE74912, GSE74912_ATACseq_All_Counts.txt.gz (590,650 peaks × 132 samples, hg19).

Sample selection (8 samples, donors as replicates):

(Bcell_3+4 and Eryth_3+4 each share a donor — minor pseudoreplication; conclusions robust; a ~ donor + condition design is the cleaner ideal.)

Prep (data-import boundary, outside Galaxy): split the matrix into 8 per-sample two-column files (peak_id <TAB> count, no header, 590,650 rows). Peak IDs are chr_start_end (e.g. chr11_60222940_60223440).

Upload: the 8 count files as a list collection ATAC counts (element identifiers exactly Bcell_1..4, Eryth_1..4 — must match the sample sheet); the sample sheet as one tabular dataset:

sample	condition	donor
Bcell_1	Bcell	Donor1022
…
Eryth_4	Erythroblast	Donor6926

4. Build sequence (all /api/tools, flat pipe-keys)

4.0 The one gotcha — DESeq2’s conditional input

DESeq2’s data selector is a conditional (select_data|how). Pass it with flat pipe-delimited keys, not a nested {select_data:{...}} dict. A nested dict mis-routes to the datasets_per_level default and silently feeds the sample sheet as the only counts file (job errors / garbage). Flat form:

select_data|how                          = sample_sheet_contrasts
select_data|countsFile                    = {src: hdca, id: <ATAC counts collection>}
select_data|sample_sheet                  = {src: hda,  id: <sample sheet>}
select_data|design_formula_mode|mode      = automatic
select_data|design_formula_mode|factor    = 2
select_data|design_formula_mode|reference_level = Erythroblast
select_data|design_formula_mode|target_level    = Bcell
header                                     = true
tximport|tximport_selector                = count
advanced_options|fit_type                 = 1
advanced_options|auto_mean_filter_off     = false
advanced_options|outlier_filter_off       = false
advanced_options|outlier_replace_off      = false
advanced_options|use_beta_priors          = false
output_options|alpha_ma                   = 0.1
output_options|output_selector            = [pdf, normCounts, normVST]

Direction: reference_level=Erythroblast ⇒ LFC > 0 = B-cell-gained, LFC < 0 = Erythroblast-gained. Confirm with a marker — MS4A1/CD20 chr11_60222940_… should be LFC ≈ +7.3.

Step 1 — DESeq2

Outputs used: result table (tabular; cols 1=peak, 2=baseMean, 3=log2FC, 4=lfcSE, 5=stat, 6=pvalue, 7=padj) and the 5-page diagnostics PDF (page 1 = PCA). Genome-wide run ≈ 6–7 min.

Step 2 — NA-filter (tp_awk) on the DESeq2 result

infile = {src: hda, id: <DESeq2 result>}
code   = $6!="NA" && $7!="NA"{print}

→ testable peaks (≈ 247,038). Prerequisite for volcano (it errors on NA-padj rows).

Step 3 — Volcano Plot on the NA-filtered table

Column params live inside the with_header conditional (top-level ones are silently ignored):

input_file                = {src: hda, id: <NA-filtered>}
with_header|header        = no
with_header|label_col     = 1
with_header|lfc_col       = 3
with_header|pval_col      = 6
with_header|fdr_col       = 7
signif_thresh = 0.05   lfc_thresh = 1.0   labels|label_select = none
plot_options|legend_labs = "Down,Not Sig,Up"

→ volcano PDF.

Step 4 — Significance filter (tp_awk) on the NA-filtered table

Branch in parallel with the volcano — do not replace the volcano’s input (it needs the full non-significant cloud).

infile = {src: hda, id: <NA-filtered>}
code   = $7+0 < 0.05 && ($3+0 >= 1 || $3+0 <= -1){print}

→ Significant differential peaks (45,620).

Step 5 — Rank (Sort) the significant peaks by log2FC, descending

infile = {src: hda, id: <significant peaks>}
header = 0
sortkeys_0|column = 3   sortkeys_0|order = r   sortkeys_0|style = g

(g = general-numeric, handles scientific notation.) Top = largest +LFC (B-cell), bottom = most −LFC (Eryth).

Step 6 — Top B-cell-gained = Select first 50 (tp_head_tool)

infile = {src: hda, id: <ranked>}   complement = ""   count = 50

Step 7 — Top Erythroblast-gained = Select last 50 (tp_tail_tool)

infile = {src: hda, id: <ranked>}   complement = ""   num_lines = 50

Rename the new outputs to meaningful names (Significant differential peaks…, Top 50 B-cell-gained peaks, Top 50 Erythroblast-gained peaks) — these become the workflow output labels.

Steps 8–11 — Figures: PDF → PNG (see §1)

Per figure (DESeq2 plots PDF, volcano PDF):

# convert
graphicsmagick_image_convert:
  input                 = {src: hda, id: <figure PDF>}
  output_format         = png
  palette|palette_select = no
  resize                = 100
→ list collection splitted_pdf (temp_000 … per page)

# extract page 1
__EXTRACT_DATASET__:
  input                 = {src: hdca, id: <splitted_pdf collection>}
  which|which_dataset   = by_identifier
  which|identifier      = temp_000
→ standalone PNG HDA

Rename the two extracted PNGs → Sample PCA plot (from DESeq2 plots page 1) and Volcano plot (from the volcano PDF).

5. Notebook page

Each directive in its own ```galaxy block, every one referencing a real on-graph output:

1. history_dataset_as_image(history_dataset_id=<Sample PCA plot PNG>) — PCA (extracted page 1 of DESeq2 plots)
2. history_dataset_as_image(history_dataset_id=<Volcano plot PNG>)
3. history_dataset_as_table(history_dataset_id=<DESeq2 result>, title="…", compact=true)
4. history_dataset_as_table(history_dataset_id=<significant peaks>, title="…", compact=true)
5. history_dataset_as_table(history_dataset_id=<top B-cell>, title="…", compact=true)
6. history_dataset_as_table(history_dataset_id=<top Eryth>, title="…", compact=true)

6. Extraction (page → workflow)

1. GET /api/pages/{page}/workflow_extraction_summary — expect 0 warnings, every tool step seeded, the 6 referenced outputs exposed.
2. POST /api/workflows/extract:

{
 "workflow_name": "Differential ATAC-seq accessibility (Eryth vs B-cell)",
 "job_ids": ["<DESeq2>","<NA-filter>","<volcano>","<sig-filter>","<sort>","<head>","<tail>","<convert PCA>","<extract PCA>","<convert volcano>","<extract volcano>"],
 "hda_ids":  ["<sample sheet>"],
 "hdca_ids": ["<ATAC counts collection>"],
 "from_page_id": "<page>",
 "report_title": "Differential ATAC-seq Accessibility — Erythroblast vs B-cell"
}

Easiest source of job_ids: the jobs[] from the summary where step_type=="tool" and checked. Pass only the collection as the counts input — NOT its 8 element HDAs (those become orphan input_dataset steps). The extractor maps DESeq2’s expanded element-inputs back to the collection.

Result: 13-step workflow, 0 dangling, 6 outputs, report markdown carried with each directive rewritten from history_dataset_id=… to workflow-relative output="<label>" (e.g. history_dataset_as_image(output="Sample PCA plot")). Each figure branch — DESeq2 → graphicsmagick_image_convert → Extract dataset and volcanoplot → graphicsmagick_image_convert → Extract dataset — is reproduced when the workflow runs.

7. Verification

8. Notes

• DESeq2 result column indices can shift with output_selector/version — confirm 3=log2FC, 6=pvalue, 7=padj from the header before wiring the awk/volcano column args.
• Figures are tool outputs (§1), so no pymupdf / core change is needed. graphicsmagick_image_convert splits a PDF into a one-page-per-element collection; temp_000 is page 1. Rename the extracted PNG before extraction so the workflow output label is meaningful (else it inherits temp_000).
• The 8 count files are prepared input data (data-import boundary), not an analysis step — legitimate workflow inputs.