feat: add germline variant calling job bundle using bioconda

job_bundles/README.md

@@ -60,6 +60,7 @@ Nested collections provide their own complete indexes where applicable.
 | [Maya V-Ray tiled render](tile_render_with_maya_vray/) | V-Ray tile rendering followed by OpenImageIO assembly | You need tiled EXR output from Maya and V-Ray |
 | [V-Ray Linux region render](tile_render_with_vray_linux/) | Region rendering, asset discovery, path mapping, and image merge | You render `.vrscene` files in parallel regions on Linux |
 | [Maya Arnold turntable](turntable_with_maya_arnold/) | Building a scene around an OBJ, then rendering frames and encoding video | You need an easy-to-submit 3D asset review utility |
+| [Germline variant calling](variant_calling_bwa/) | DNA sequence alignment and variant calling with bioconda software, scattered per sample and per genome region | You want to find genetic differences from sequencing data, or a sequence-bioinformatics starting point |
 | [Studio VFX shot pipeline](vfx_pipeline/) | A launcher, static render bundle, Conda software delivery, and Flow publishing mapped onto a studio's existing pipeline | You are porting a studio shot pipeline to Deadline Cloud |
 | [AutoDock Vina virtual screening](virtual_screening_vina/) | Parallel molecular docking and ranked result aggregation | You want to screen many ligands against a protein target |
 | [vLLM evaluation leaderboard](vllm_lm_eval_leaderboard/) | Parallel model evaluation and final CSV/Markdown aggregation | You need to compare multiple LLMs across benchmarks |

job_bundles/variant_calling_bwa/sample_inputs/.gitignore

@@ -0,0 +1,5 @@
+# Test data is downloaded by fetch_test_data.py, not committed. The repository
+# root .gitignore covers *.tar.gz and *.zip archives but not the bare *.fastq.gz
+# and *.fasta files fetched here, so exclude both directories outright.
+reads/
+reference/

job_bundles/variant_calling_bwa/sample_inputs/README.md

@@ -0,0 +1,41 @@
+# Sample inputs
+
+The data this sample runs on is downloaded rather than committed, because test data is better
+fetched from its upstream source than vendored. A `.gitignore` in this directory keeps the
+downloads out of version control.
+
+```console
+python fetch_test_data.py
+```
+
+That writes about 950 KB:
+
+```
+reads/
+  tiny_n_R1.fastq.gz   tiny_n_R2.fastq.gz    # "normal" library
+  tiny_t_R1.fastq.gz   tiny_t_R2.fastq.gz    # "tumor" library
+reference/
+  human_g1k_v37_decoy.small.fasta            # 6-contig GRCh37 subset
+  human_g1k_v37_decoy.small.fasta.fai
+```
+
+The read sets are the normal and tumor libraries of a synthetic pair, but because this sample does no
+somatic calling, they act as two independent samples to fan out over.
+
+The reference uses GRCh37-style contig naming with no `chr` prefix, and `fetch_test_data.py` prints the
+contig names once the download finishes. The reads align to only part of contig `1`, so the job
+template's default `Regions` is a set of windows of that contig rather than the contig names.
+
+Only lane `L001` of each library is downloaded. The upstream location also has `L002`, if you want
+roughly twice the reads.
+
+## Provenance
+
+The reads and reference are provided by
+[AWS HealthOmics for their tutorials](https://github.com/aws-samples/aws-healthomics-tutorials),
+hosted in the public, unauthenticated `aws-genomics-static-us-east-1` bucket under
+`omics-data/test-datasets/nf-core-sarek/`. The data originates from
+[nf-core/test-datasets](https://github.com/nf-core/test-datasets) and is MIT licensed. A copy of
+that license is in the bucket alongside the data.
+
+`fetch_test_data.py --clean` removes what it previously downloaded before fetching again.

job_bundles/variant_calling_bwa/sample_inputs/fetch_test_data.py

@@ -0,0 +1,142 @@
+#!/usr/bin/env python3
+"""Download the tiny public test dataset used by this sample.
+
+The reads and reference are provided by AWS HealthOmics for their tutorials
+(https://github.com/aws-samples/aws-healthomics-tutorials), hosted in the public,
+unauthenticated ``aws-genomics-static-us-east-1`` bucket. The data originates
+from nf-core/test-datasets and is MIT licensed. The whole set is under 1 MB, so
+the pipeline runs end to end in seconds.
+
+Reads land in ``reads/`` renamed to the ``<sample>_R1.fastq.gz`` /
+``<sample>_R2.fastq.gz`` convention the job template expects. The reference and
+its .fai land in ``reference/``.
+
+Usage:
+    python fetch_test_data.py            # download into ./reads and ./reference
+    python fetch_test_data.py --clean    # remove downloaded files first
+"""
+from __future__ import annotations
+
+import argparse
+import shutil
+import sys
+import urllib.error
+import urllib.request
+from pathlib import Path
+
+BASE_URL = (
+    "https://aws-genomics-static-us-east-1.s3.amazonaws.com"
+    "/omics-data/test-datasets/nf-core-sarek"
+)
+
+# One lane per sample keeps the demo small. tiny_n is the "normal" library and
+# tiny_t the "tumor" library of the same synthetic pair; here they simply act as
+# two independent samples to fan out over.
+READS = {
+    "tiny_n_R1.fastq.gz": "testdata/tiny/normal/tiny_n_L001_R1_xxx.fastq.gz",
+    "tiny_n_R2.fastq.gz": "testdata/tiny/normal/tiny_n_L001_R2_xxx.fastq.gz",
+    "tiny_t_R1.fastq.gz": "testdata/tiny/tumor/tiny_t_L001_R1_xxx.fastq.gz",
+    "tiny_t_R2.fastq.gz": "testdata/tiny/tumor/tiny_t_L001_R2_xxx.fastq.gz",
+}
+
+# The .fai is fetched so BuildIndex can skip samtools faidx. The bwa index files
+# are deliberately not fetched: BuildIndex rebuilds them in a second or two for a
+# reference this small, which keeps the download to a single FASTA and avoids
+# depending on the index having been built by a compatible bwa version.
+REFERENCE = {
+    "human_g1k_v37_decoy.small.fasta": "reference/human_g1k_v37_decoy.small.fasta",
+    "human_g1k_v37_decoy.small.fasta.fai": "reference/human_g1k_v37_decoy.small.fasta.fai",
+}
+
+HERE = Path(__file__).resolve().parent
+
+
+def download(url: str, dest: Path) -> None:
+    if dest.exists() and dest.stat().st_size > 0:
+        print(f"  exists, skipping: {dest.name}")
+        return
+    dest.parent.mkdir(parents=True, exist_ok=True)
+    print(f"  {dest.name} <- {url}")
+
+    # Download to a temporary name and rename only once the transfer is complete
+    # and the byte count matches Content-Length. A partial file left at the final
+    # path would be skipped by the check above on every later run, and the failure
+    # would not surface until a tool choked on the truncated data.
+    partial = dest.with_name(dest.name + ".part")
+    try:
+        with urllib.request.urlopen(url, timeout=120) as response:
+            expected = response.headers.get("Content-Length")
+            with partial.open("wb") as out:
+                shutil.copyfileobj(response, out)
+        written = partial.stat().st_size
+        if expected is not None and written != int(expected):
+            raise OSError(f"expected {int(expected):,} bytes, received {written:,}")
+        partial.replace(dest)
+    except urllib.error.HTTPError as exc:
+        partial.unlink(missing_ok=True)
+        sys.exit(f"ERROR: HTTP {exc.code} fetching {url}")
+    except urllib.error.URLError as exc:
+        partial.unlink(missing_ok=True)
+        sys.exit(f"ERROR: could not reach {url}: {exc.reason}")
+    except (OSError, KeyboardInterrupt):
+        partial.unlink(missing_ok=True)
+        raise
+    print(f"      {dest.stat().st_size:,} bytes")
+
+
+def main() -> int:
+    parser = argparse.ArgumentParser(
+        description=__doc__,
+        formatter_class=argparse.RawDescriptionHelpFormatter,
+    )
+    parser.add_argument(
+        "--clean",
+        action="store_true",
+        help="delete the reads/ and reference/ contents before downloading",
+    )
+    args = parser.parse_args()
+
+    reads_dir = HERE / "reads"
+    reference_dir = HERE / "reference"
+
+    if args.clean:
+        for directory in (reads_dir, reference_dir):
+            if directory.exists():
+                print(f"Removing {directory}")
+                shutil.rmtree(directory)
+
+    print("Downloading reads (nf-core/test-datasets, MIT licensed):")
+    for name, key in READS.items():
+        download(f"{BASE_URL}/{key}", reads_dir / name)
+
+    print("Downloading reference:")
+    for name, key in REFERENCE.items():
+        download(f"{BASE_URL}/{key}", reference_dir / name)
+
+    # Count only the downloaded data, not this script and the docs beside it.
+    total = sum(
+        f.stat().st_size
+        for directory in (reads_dir, reference_dir)
+        for f in directory.rglob("*")
+        if f.is_file()
+    )
+    print(f"\nDone. {total / 1024:.0f} KB downloaded into {HERE}")
+
+    fai = reference_dir / "human_g1k_v37_decoy.small.fasta.fai"
+    if fai.exists():
+        names = [
+            line.split("\t")[0]
+            for line in fai.read_text(encoding="utf-8").splitlines()
+            if line.strip()
+        ]
+        print(f"Reference contigs: {','.join(names)}")
+        print(
+            f"The sample reads align to only part of contig {names[0]}, which is "
+            "why the job template's default Regions is a set of windows within it "
+            "rather than these contig names."
+        )
+    return 0
+
+
+if __name__ == "__main__":
+    sys.exit(main())

job_bundles/variant_calling_bwa/scripts/align_sample.sh

@@ -0,0 +1,86 @@
+#!/usr/bin/env bash
+# Align one sample with bwa mem, then sort and index the BAM.
+set -euo pipefail
+source "$(dirname "${BASH_SOURCE[0]}")/common.sh"
+
+SAMPLES_CSV=""
+SAMPLE_INDEX=""
+READS=""
+OUT=""
+REFERENCE_NAME=""
+WORK_ROOT=""
+for arg in "$@"; do
+  case "$arg" in
+    --samples=*)         SAMPLES_CSV="${arg#*=}" ;;
+    --sample-index=*)    SAMPLE_INDEX="${arg#*=}" ;;
+    --reads-dir=*)       READS="${arg#*=}" ;;
+    --output-dir=*)      OUT="${arg#*=}" ;;
+    --reference-name=*)  REFERENCE_NAME="${arg#*=}" ;;
+    --session-work-dir=*) WORK_ROOT="${arg#*=}" ;;
+    *) unknown_arg "$arg" ;;
+  esac
+done
+require_arg "$SAMPLES_CSV" --samples
+require_arg "$SAMPLE_INDEX" --sample-index
+require_arg "$READS" --reads-dir
+require_arg "$OUT" --output-dir
+require_arg "$REFERENCE_NAME" --reference-name
+require_arg "$WORK_ROOT" --session-work-dir
+
+parse_list SAMPLES "$SAMPLES_CSV"
+require_non_empty_list "${#SAMPLES[@]}" Samples
+validate_sample_names "${SAMPLES[@]}"
+if (( SAMPLE_INDEX >= ${#SAMPLES[@]} )); then
+  echo "Sample index $SAMPLE_INDEX is past the end of the sample list; nothing to do."
+  exit 0
+fi
+SAMPLE="${SAMPLES[$SAMPLE_INDEX]}"
+BAM_DIR="$OUT/alignments"
+# Named by index rather than sample name: this directory is removed with 'rm -rf'
+# below, so keeping the name a plain integer means the target cannot depend on
+# the contents of a job parameter.
+WORK="$WORK_ROOT/align_${SAMPLE_INDEX}"
+mkdir -p "$BAM_DIR" "$WORK"
+
+R1="$READS/${SAMPLE}_R1.fastq.gz"
+R2="$READS/${SAMPLE}_R2.fastq.gz"
+BAM="$BAM_DIR/${SAMPLE}.sorted.bam"
+
+# Derived from the output directory, not read from a file: BuildIndex ran in a
+# different session directory.
+IDX_PREFIX="$OUT/reference/$REFERENCE_NAME"
+# bwa mem takes the prefix and finds the index itself, accepting either the
+# "<reference>.bwt" or the "<reference>.64.bwt" naming, so accept both here too.
+if [[ ! -f "${IDX_PREFIX}.bwt" && ! -f "${IDX_PREFIX}.64.bwt" ]]; then
+  echo "ERROR: no bwa index beside ${IDX_PREFIX}. Did BuildIndex run?" >&2
+  exit 1
+fi
+echo "Aligning $SAMPLE against $IDX_PREFIX"
+
+# Scale to the worker rather than hardcoding a thread count. The step requires a
+# minimum of 2 vCPU, not exactly 2, so a fleet is free to place this task on a
+# much larger instance; a fixed count would leave most of it idle.
+ALIGN_THREADS="$(nproc)"
+# 'samtools sort -m' is per thread, so the sort's total appetite is threads times
+# -m. Divide a fixed budget between the threads to keep that total flat, and cap
+# the thread count so it stays flat on a very large worker: sorting is I/O bound
+# well before this many threads anyway.
+SORT_THREADS=$(( ALIGN_THREADS > 1 ? ALIGN_THREADS / 2 : 1 ))
+(( SORT_THREADS > 4 )) && SORT_THREADS=4
+SORT_MEM_MB=$(( 4096 / SORT_THREADS ))
+echo "Using $ALIGN_THREADS alignment thread(s), $SORT_THREADS sort thread(s) at ${SORT_MEM_MB}M each"
+
+# The @RG line carries the sample name into the BAM so that bcftools attributes
+# calls to the right sample column in the VCF.
+bwa mem \
+  -t "$ALIGN_THREADS" \
+  -R "@RG\tID:${SAMPLE}\tSM:${SAMPLE}\tPL:ILLUMINA\tLB:${SAMPLE}" \
+  "$IDX_PREFIX" "$R1" "$R2" \
+  | samtools sort -@ "$SORT_THREADS" -m "${SORT_MEM_MB}M" -T "$WORK/sort" -o "$BAM" -
+
+samtools index "$BAM"
+samtools flagstat "$BAM" > "$BAM_DIR/${SAMPLE}.flagstat.txt"
+echo "--- flagstat: $SAMPLE ---"
+head -n 5 "$BAM_DIR/${SAMPLE}.flagstat.txt"
+rm -rf "$WORK"
+echo "Alignment complete: $BAM"

job_bundles/variant_calling_bwa/scripts/build_index.sh

@@ -0,0 +1,100 @@
+#!/usr/bin/env bash
+# Create the .fai and bwa indexes once, if they are not already present.
+set -euo pipefail
+source "$(dirname "${BASH_SOURCE[0]}")/common.sh"
+
+REFERENCE_DIR=""
+REFERENCE_NAME=""
+OUT=""
+for arg in "$@"; do
+  case "$arg" in
+    --reference-dir=*)  REFERENCE_DIR="${arg#*=}" ;;
+    --reference-name=*) REFERENCE_NAME="${arg#*=}" ;;
+    --output-dir=*)     OUT="${arg#*=}" ;;
+    *) unknown_arg "$arg" ;;
+  esac
+done
+require_arg "$REFERENCE_DIR" --reference-dir
+require_arg "$REFERENCE_NAME" --reference-name
+require_arg "$OUT" --output-dir
+
+REF="$REFERENCE_DIR/$REFERENCE_NAME"
+mkdir -p "$OUT"
+
+# The reference name is a name within the reference directory, so a value
+# containing a separator would resolve outside the staged directory.
+if [[ "$REFERENCE_NAME" != "$(basename "$REFERENCE_NAME")" ]]; then
+  echo "ERROR: ReferenceFastaName must be a filename, not a path: $REFERENCE_NAME" >&2
+  exit 1
+fi
+
+if [[ ! -f "$REF" ]]; then
+  echo "ERROR: reference FASTA not found: $REF" >&2
+  echo "ReferenceDir contains:" >&2
+  ls -1 "$REFERENCE_DIR" >&2 || true
+  exit 1
+fi
+
+# Assemble the reference and both indexes under $OUT/reference. Later steps
+# rebuild this location from the OutputDir parameter rather than reading a path
+# from here, because each step runs in its own session directory. Copying also
+# leaves the staged input untouched, which matters when job attachments stages
+# it read-only.
+REF_DIR="$OUT/reference"
+LOCAL_REF="$REF_DIR/$REFERENCE_NAME"
+mkdir -p "$REF_DIR"
+cp -f "$REF" "$LOCAL_REF"
+
+if [[ -f "${REF}.fai" ]]; then
+  cp -f "${REF}.fai" "${LOCAL_REF}.fai"
+  echo "Reused the .fai staged beside the reference."
+else
+  samtools faidx "$LOCAL_REF"
+  echo "Built ${LOCAL_REF}.fai"
+fi
+
+# A bwa index is a set of files sharing the reference's prefix, and bwa mem needs
+# all of them, so a partial set is rebuilt rather than copied.
+#
+# 'bwa index -6' names them "<reference>.64.*" instead of "<reference>.*", and bwa
+# looks for that variant first, so check for it the same way.
+BWA_EXTS=(amb ann bwt pac sa)
+
+# Echoes the infix of a complete index beside $1, or nothing if neither is complete.
+complete_bwa_index_infix() {
+  local _prefix="$1" _infix _ext _complete
+  for _infix in ".64" ""; do
+    _complete=1
+    for _ext in "${BWA_EXTS[@]}"; do
+      [[ -f "${_prefix}${_infix}.${_ext}" ]] || _complete=0
+    done
+    if [[ "$_complete" -eq 1 ]]; then
+      printf '%s' "$_infix"
+      return 0
+    fi
+  done
+  return 1
+}
+
+# Clear both variants first. Reusing an OutputDir can leave an index from an
+# earlier run here, and bwa looks for the ".64" one before the standard one, so a
+# leftover ".64" set would be loaded in preference to whichever index this run
+# just put in place.
+for infix in ".64" ""; do
+  for ext in "${BWA_EXTS[@]}"; do
+    rm -f "${LOCAL_REF}${infix}.${ext}"
+  done
+done
+
+if BWA_INFIX="$(complete_bwa_index_infix "$REF")"; then
+  for ext in "${BWA_EXTS[@]}"; do
+    cp -f "${REF}${BWA_INFIX}.${ext}" "${LOCAL_REF}${BWA_INFIX}.${ext}"
+  done
+  echo "Reused the bwa index staged beside the reference."
+else
+  bwa index "$LOCAL_REF"
+  echo "Built the bwa index for $LOCAL_REF"
+fi
+
+echo "Reference ready at $LOCAL_REF"
+ls -1 "$REF_DIR"