Skip to content

Paired Experiment-Control (Multiple Pairs)#

The scalable version of gallery 14: the same somatic pipeline driven by [[pairs]] with {experiment} / {control} / {pair_id} wildcards, fanning out one branch per pair.

What It Demonstrates#

  • [[pairs]] entries expand wildcard rules per pair; aggregation rules without pair wildcards stay single-instance
  • --germline-resource / --panel-of-normals are noted in-rule for production hardening
  • Read-group SM tags per branch make Mutect2 -normal {control} resolve correctly

Workflow Definition#

# examples/gallery/15_paired_experiment_control_pairs.oxoflow
# 15 — Paired Experiment-Control (Multiple Pairs)
# The scalable version of 14: the same somatic pipeline driven by
# [[pairs]] with {experiment}/{control}/{pair_id} wildcards.

[workflow]
name = "multi-case-experiment-control"
version = "1.0.0"
description = "Parallel somatic variant calling for multiple experiment-control pairs"
author = "oxo-flow examples"

[config]
reference_fasta = "/data/references/GRCh38/genome.fa"
genome_build = "GRCh38"
known_sites = "/data/references/GRCh38/dbsnp_146.hg38.vcf.gz"

[defaults]
threads = 8
memory = "16G"

[[rules]]
name = "fastp_experiment"
input = [
    "raw/{experiment}_R1.fq.gz",
    "raw/{experiment}_R2.fq.gz",
]
output = [
    "trimmed/{experiment}_R1.fq.gz",
    "trimmed/{experiment}_R2.fq.gz",
    "qc/{experiment}_fastp.json",
]
shell = "fastp -i {input[0]} -I {input[1]} -o {output[0]} -O {output[1]} --json {output[2]} --thread {threads}"

[rules.resources]
threads = 8

[rules.environment]
conda = "envs/fastp.yaml"

[[rules]]
name = "fastp_control"
input = [
    "raw/{control}_R1.fq.gz",
    "raw/{control}_R2.fq.gz",
]
output = [
    "trimmed/{control}_R1.fq.gz",
    "trimmed/{control}_R2.fq.gz",
    "qc/{control}_fastp.json",
]
shell = "fastp -i {input[0]} -I {input[1]} -o {output[0]} -O {output[1]} --json {output[2]} --thread {threads}"

[rules.resources]
threads = 8

[rules.environment]
conda = "envs/fastp.yaml"

[[rules]]
name = "bwa_mem2_experiment"
input = [
    "trimmed/{experiment}_R1.fq.gz",
    "trimmed/{experiment}_R2.fq.gz",
]
output = ["aligned/{experiment}.sorted.bam"]
shell = """
rg=$(printf '@RG\tID:%s\tSM:%s\tPL:ILLUMINA' '{experiment}' '{experiment}')
bwa-mem2 mem -t {threads} -R "$rg" {config.reference_fasta} {input[0]} {input[1]} | samtools sort -@ 4 -o {output[0]}
"""

[rules.resources]
threads = 16
memory = "32G"

[rules.environment]
conda = "envs/bwa_mem2.yaml"

[[rules]]
name = "bwa_mem2_control"
input = [
    "trimmed/{control}_R1.fq.gz",
    "trimmed/{control}_R2.fq.gz",
]
output = ["aligned/{control}.sorted.bam"]
shell = """
rg=$(printf '@RG\tID:%s\tSM:%s\tPL:ILLUMINA' '{control}' '{control}')
bwa-mem2 mem -t {threads} -R "$rg" {config.reference_fasta} {input[0]} {input[1]} | samtools sort -@ 4 -o {output[0]}
"""

[rules.resources]
threads = 16
memory = "32G"

[rules.environment]
conda = "envs/bwa_mem2.yaml"

[[rules]]
name = "markdup_experiment"
input = ["aligned/{experiment}.sorted.bam"]
output = [
    "dedup/{experiment}.dedup.bam",
    "dedup/{experiment}.metrics.txt",
]
shell = "gatk MarkDuplicates --VALIDATION_STRINGENCY SILENT -I {input[0]} -O {output[0]} -M {output[1]}"

[rules.resources]
threads = 4
memory = "16G"

[rules.environment]
conda = "envs/gatk.yaml"

[[rules]]
name = "markdup_control"
input = ["aligned/{control}.sorted.bam"]
output = [
    "dedup/{control}.dedup.bam",
    "dedup/{control}.metrics.txt",
]
shell = "gatk MarkDuplicates --VALIDATION_STRINGENCY SILENT -I {input[0]} -O {output[0]} -M {output[1]}"

[rules.resources]
threads = 4
memory = "16G"

[rules.environment]
conda = "envs/gatk.yaml"

[[rules]]
name = "bqsr_experiment"
input = ["dedup/{experiment}.dedup.bam"]
output = ["recal/{experiment}.recal.bam"]
shell = """
  gatk BaseRecalibrator --read-validation-stringency SILENT -I {input[0]} -R {config.reference_fasta} --known-sites {config.known_sites} -O recal/{experiment}.recal.table && gatk ApplyBQSR -I {input[0]} -R {config.reference_fasta} --bqsr-recal-file recal/{experiment}.recal.table -O {output[0]}
"""

[rules.resources]
threads = 4
memory = "16G"

[rules.environment]
conda = "envs/gatk.yaml"

[[rules]]
name = "bqsr_control"
input = ["dedup/{control}.dedup.bam"]
output = ["recal/{control}.recal.bam"]
shell = """
  gatk BaseRecalibrator --read-validation-stringency SILENT -I {input[0]} -R {config.reference_fasta} --known-sites {config.known_sites} -O recal/{control}.recal.table && gatk ApplyBQSR -I {input[0]} -R {config.reference_fasta} --bqsr-recal-file recal/{control}.recal.table -O {output[0]}
"""

[rules.resources]
threads = 4
memory = "16G"

[rules.environment]
conda = "envs/gatk.yaml"

[[rules]]
name = "mutect2"
input = [
    "recal/{experiment}.recal.bam",
    "recal/{control}.recal.bam",
]
output = ["variants/{pair_id}.mutect2.vcf.gz"]
shell = """
  # Production runs should also pass --germline-resource and --panel-of-normals
  gatk Mutect2 --read-validation-stringency SILENT -I {input[0]} -I {input[1]} -normal {control} -R {config.reference_fasta} -O {output[0]}
"""

[rules.resources]
threads = 4
memory = "16G"

[rules.environment]
conda = "envs/gatk.yaml"

[[rules]]
name = "filter_mutect_calls"
input = ["variants/{pair_id}.mutect2.vcf.gz"]
output = ["variants/{pair_id}.mutect2.filtered.vcf.gz"]
shell = "gatk FilterMutectCalls -V {input[0]} -R {config.reference_fasta} -O {output[0]}"

[rules.resources]
threads = 2
memory = "8G"

[rules.environment]
conda = "envs/gatk.yaml"

[[rules]]
name = "annotate_somatic"
input = ["variants/{pair_id}.mutect2.filtered.vcf.gz"]
output = ["annotated/{pair_id}.annotated.vcf.gz"]
shell = "vep --input_file {input[0]} --output_file {output[0]} --format vcf --vcf --offline --cache"

[rules.resources]
threads = 4
memory = "8G"

[rules.environment]
conda = "envs/vep.yaml"

[[rules]]
name = "clinical_report"
input = ["annotated/{pair_id}.annotated.vcf.gz"]
output = ["reports/{pair_id}_clinical_report.html"]
shell = "python scripts/generate_report.py --input {input[0]} --output {output[0]} --pair {pair_id}"
description = "Generate clinical-grade variant report for a experiment-control pair"

[rules.resources]
threads = 2
memory = "4G"

[rules.environment]
conda = "envs/report.yaml"

[[pairs]]
pair_id = "CASE_001"
experiment = "EXP_01"
control = "CTRL_01"
experiment_type = "lung_adenocarcinoma"

[pairs.metadata]

[[pairs]]
pair_id = "CASE_002"
experiment = "EXP_02"
control = "CTRL_02"
experiment_type = "colorectal"

[pairs.metadata]

Try It#

# Inspect the expanded plan first — no data needed:
oxo-flow dry-run examples/gallery/15_paired_experiment_control_pairs.oxoflow

# Copy the environment specs next to the workflow, adapt [config]
# paths to your data, then run:
oxo-flow run examples/gallery/15_paired_experiment_control_pairs.oxoflow

Input data and environments

Input paths under /data/references/... and raw/ are placeholders — replace them with your own data. The referenced envs/*.yaml specs ship in examples/envs/.