diff --git a/job_bundles/README.md b/job_bundles/README.md index e589ff39..1d972230 100644 --- a/job_bundles/README.md +++ b/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 | diff --git a/job_bundles/variant_calling_bwa/README.md b/job_bundles/variant_calling_bwa/README.md new file mode 100644 index 00000000..4ebedb51 --- /dev/null +++ b/job_bundles/variant_calling_bwa/README.md @@ -0,0 +1,494 @@ +# Germline variant calling with bwa, samtools, and bcftools + +Finds the genetic differences between an individual and a reference genome. + +A DNA sequencing machine does not read a genome end to end. It produces hundreds of millions of +short fragments, each a few hundred letters long, from random positions. Turning those fragments +into a list of differences takes three stages, and this job bundle runs all three: + +1. **Align.** Work out where each fragment came from by matching it against a reference genome. + `bwa` does this. +2. **Call.** At every position, compare the stacked-up fragments to the reference and decide whether + the difference is real or a sequencing error. `bcftools` does this. +3. **Merge.** Collect the results into one file listing every confident difference. + +The differences are **variants**, and finding them is **variant calling**, the workhorse analysis +behind everything from diagnosing inherited disease to breeding drought-tolerant crops. *Germline* +means the variants an individual was born with, rather than ones a tumor acquired later. The input +format is **FASTQ** (the sequencer's fragments) and the output a **VCF** (the list of variants). + +Both expensive stages divide cleanly, which is what suits this to Deadline Cloud. Each sample aligns +on its own worker, and separate stretches of the genome are then called independently and stitched +back together. + +## What this sample demonstrates + +* **Software from bioconda.** `bwa`, `samtools`, `bcftools`, and `fastqc` are declared as conda + packages rather than built into container images. See + [Conda channel order](#conda-channel-order) for why `conda-forge` must come first. +* **Fan-out over samples and over genome regions.** `AlignReads` fans out over samples and + `CallVariants` over regions, using different task parameters in different steps. `CallVariants` + is the scatter and `MergeVariants` the gather. +* **A job environment as a precondition check.** `BioToolchain` verifies every command in + `RequiredTools` is on `PATH` once per session and, if any is missing, prints the exact packages and + channels to configure. +* **Scripts bundled beside the template.** The shell lives in [`scripts/`](scripts/) rather than + embedded in `template.yaml`, following stage 2 of the + [job development progression](../job_dev_progression/). See + [Bundled scripts](#bundled-scripts) for what that changes. + +The tool sequence comes from the +[AWS HealthOmics WDL variant-calling tutorial pipeline](https://github.com/aws-samples/aws-healthomics-tutorials/tree/main/example-workflows/wdl/variant-calling-pipeline), +reimplemented here for Open Job Description. [WDL](https://openwdl.org/) and Open Job Description are +both open workflow specifications with independent implementations, so having the same pipeline in +both is a good way to compare their designs. WDL infers the order of work from how data +flows between tasks and delivers software as containers; this bundle declares step dependencies +explicitly and installs its tools from conda. Either way the tool builds are equivalent, since the +BioContainers images WDL uses are themselves built from bioconda packages. + +This pipeline also differs from the original in two ways. That one scatters over samples and gathers +them into a single whole-genome call; this one adds a second scatter over regions, so calling is +distributed rather than left on one worker. It also takes paired-end reads, the form most sequencers +produce, where the original takes one FASTQ per sample. + +## How it works + +``` +QualityControl (1 task/sample) BuildIndex (1 task) + fastqc on each read pair samtools faidx + bwa index + → output/qc/ → output/reference/ + no dependencies, starts at once │ + │ ▼ + │ AlignReads (1 task/sample) + │ bwa mem | samtools sort, then index + │ → output/alignments/.sorted.bam + │ │ + │ ▼ + │ CallVariants (1 task/region) ── the scatter + │ bcftools mpileup | call | filter, + │ all samples jointly per region + │ → output/vcf_by_region/ + │ │ + └────────────────┬─────────────────┘ + ▼ + MergeVariants (1 task) ── the gather + bcftools concat -d exact | norm, stats, MultiQC + → output/variants.vcf.gz +``` + +The default parameters (2 samples, 4 regions) produce 10 tasks. + +The default regions are four windows of contig `1` rather than whole contigs, because the sample reads +align to only part of it and scattering over every contig would leave most tasks with nothing to call. +Windowing is also how you scatter a real genome: chromosomes vary widely in size, so a per-chromosome +scatter leaves one task running long after the rest finish. Keep windows non-overlapping, since +overlapping ones call the same site twice. + +`RegionRange` must cover every entry in `Regions`. The practical ceiling is about 100 entries, because +a job parameter string is capped at 1024 characters. Staging the list as a data file instead lifts that +ceiling; see [Scaling the region list past the parameter limit](#scaling-the-region-list-past-the-parameter-limit). + +`CallVariants` calls every sample jointly per region. Joint calling lets the caller distinguish a site +that matches the reference in one sample from one that merely lacked coverage there. + +## What this sample leaves out + +A production analysis would add steps this bundle omits to stay readable: + +* **Duplicate marking** (`samtools fixmate -m` then `samtools markdup`, or Picard). PCR and optical + duplicates inflate apparent allele support. `bcftools mpileup` already skips reads flagged as + duplicates, but nothing here sets that flag, so `flagstat` reports `0 duplicates`. The + `LB:` field in the read group exists for this step. +* **Base quality score recalibration.** A GATK idiom with no direct bcftools equivalent. +* **Adapter and quality trimming.** `bwa mem` soft-clips adapters, so this matters less than it once + did. Note that `QualityControl` runs alongside `AlignReads` rather than before it, so the FastQC + report cannot gate alignment. +* **Per-contig ploidy.** `Ploidy` applies job-wide. Real analyses need `bcftools call --ploidy-file` + so chrX, chrY, and the mitochondrion are treated correctly. + +Filtering and normalization are the exception, included because a callset is unusable without them: +`CallVariants` applies `MinQual` and `MinDepth`, and `MergeVariants` runs `bcftools norm` to +left-align indels and split multiallelic records. + +One scaling limit worth knowing: `bcftools mpileup` joint calling has no gVCF equivalent, so adding a +sample to a cohort means recalling all of them. + +## Prerequisites + +1. **A Deadline Cloud farm with a Linux fleet.** The tools are only published for Linux and + macOS, so the steps declare `attr.worker.os.family: linux`. + +2. **A conda queue environment** whose channels include `conda-forge` and `bioconda`. See the + [queue environment samples](../../queue_environments/) and + [Create a queue environment](https://docs.aws.amazon.com/deadline-cloud/latest/userguide/create-queue-environment.html). + If the queue is missing the tools, the job fails on entry with the packages and channels it + needs rather than partway through a step. + +3. **The Deadline Cloud CLI:** + ```console + pip install deadline + ``` + +### Conda channel order + +Set the channels to `conda-forge bioconda`, in that order. bioconda packages depend on +conda-forge for their runtime libraries, and [bioconda's documentation](https://bioconda.github.io/) +specifies conda-forge at higher priority, with strict channel priority. Reversing the order causes +dependency resolution failures that can be hard to read. + +## Setup + +Every command below runs from this sample's directory: + +```console +cd job_bundles/variant_calling_bwa +``` + +The sample data is not committed to this repository. Download it with: + +```console +python sample_inputs/fetch_test_data.py +``` + +That fetches about 950 KB: two paired-end read sets and a small GRCh37 subset as the reference. The +data is provided by +[AWS HealthOmics for their tutorials](https://github.com/aws-samples/aws-healthomics-tutorials) in the +public `aws-genomics-static-us-east-1` bucket, and originates from +[nf-core/test-datasets](https://github.com/nf-core/test-datasets), MIT licensed. See +[`sample_inputs/README.md`](sample_inputs/README.md) for what lands where. + +Bring your own data instead by pointing `ReadsDir` at a directory of +`_R1.fastq.gz` / `_R2.fastq.gz` pairs, `ReferenceDir` at the directory holding your +reference, and `ReferenceFastaName` at the FASTA filename within it, then setting `Samples` and +`Regions` to match. Remember to update `SampleRange` and `RegionRange` as described below. + +If that directory already holds a `.fai` or a bwa index built beside the FASTA, `BuildIndex` reuses +them instead of rebuilding, which is worth minutes on a whole genome. The whole bwa index must be +there to be reused, in either the `.*` or the `.64.*` naming that +`bwa index -6` produces. An incomplete set is rebuilt rather than copied. + +## Run or submit + +```console +# Submit with the defaults +deadline bundle submit . + +# Review parameters in a GUI first +deadline bundle gui-submit . + +# A fast smoke test: one sample, one region covering all the sample reads +deadline bundle submit . \ + -p Samples=tiny_n -p SampleRange=0 \ + -p Regions=1:131000-140999 -p RegionRange=0 + +# Your own cohort across whole chromosomes +deadline bundle submit . \ + -p ReadsDir=/data/cohort/reads \ + -p ReferenceDir=/data/refs \ + -p ReferenceFastaName=Homo_sapiens_assembly38.fasta \ + -p Samples=NA12878,NA12891,NA12892 -p SampleRange=0-2 \ + -p Regions=chr20,chr21,chr22 -p RegionRange=0-2 +``` + +Download the results once the job finishes: + +```console +deadline job download-output --job-id +``` + +### Run it locally with the Open Job Description CLI + +The [Open Job Description CLI](https://github.com/OpenJobDescription/openjd-cli) runs the template +without a farm, which is the fastest way to iterate. On a Linux or macOS host with the tools +available: + +```console +# Validate the template and inspect the task graph +openjd check template.yaml +openjd summary template.yaml + +# Run one sample and one region end to end, including the steps it depends on +openjd run template.yaml --step CallVariants --run-dependencies \ + -p Samples=tiny_n -p SampleRange=0 \ + -p Regions=1:131000-140999 -p RegionRange=0 + +# Run everything +openjd run template.yaml +``` + +Create the environment the tools come from with: + +```console +conda create -n variant-calling -c conda-forge -c bioconda \ + bwa samtools bcftools fastqc multiqc +conda activate variant-calling +``` + +## Parameters and outputs + +| Parameter | Default | Description | +|---|---|---| +| `ReadsDir` | `sample_inputs/reads` | Directory of `_R1.fastq.gz` / `_R2.fastq.gz` pairs | +| `ReferenceDir` | `sample_inputs/reference` | Directory holding the reference FASTA and any indexes beside it | +| `ReferenceFastaName` | `human_g1k_v37_decoy.small.fasta` | FASTA filename within the reference directory | +| `OutputDir` | `output` | Destination for all results | +| `Samples` | `tiny_n,tiny_t` | Comma-separated sample names to align in parallel | +| `SampleRange` | `0-1` | Which `Samples` entries to align, as indices from 0 | +| `Regions` | four windows of contig `1` | Comma-separated regions to call in parallel | +| `RegionRange` | `0-3` | Which `Regions` entries to call, as indices from 0 | +| `MinMappingQuality` | `20` | `bcftools mpileup -q` | +| `MinBaseQuality` | `20` | `bcftools mpileup -Q` | +| `Ploidy` | `2` | Ploidy preset for `bcftools call`; job-wide | +| `MinQual` | `20.0` | Discard called sites below this QUAL; 0 keeps everything | +| `MinDepth` | `5` | Discard called sites below this `INFO/DP`; 0 keeps everything | +| `CondaPackages` | `bwa samtools bcftools fastqc multiqc` | Packages the queue environment installs | +| `CondaChannels` | `conda-forge bioconda` | Channels, highest priority first | +| `RequiredTools` | `bwa,samtools,bcftools,fastqc,multiqc` | Commands `BioToolchain` requires on `PATH` | + +Outputs, all under `OutputDir`: + +| Path | Contents | +|---|---| +| `variants.vcf.gz` (+ `.tbi`) | The merged, deduplicated, normalized VCF, the main result | +| `variant_summary.txt` | Variant counts per contig and `bcftools stats` output | +| `vcf_by_region/region_NNNN_.vcf.gz` | Per-region calls, one per scatter task | +| `alignments/.sorted.bam` (+ `.bai`) | Sorted alignments per sample | +| `alignments/.flagstat.txt` | Alignment summary per sample | +| `qc/` | FastQC reports per read file | +| `multiqc/multiqc_report.html` | Aggregate QC report | +| `reference/` | The reference plus its `.fai` and bwa indexes, assembled by `BuildIndex` | + +### Bundled scripts + +The shell for every step lives in [`scripts/`](scripts/), not inline in `template.yaml`, following +stage 2 of the [job development progression](../job_dev_progression/). At five steps and roughly +400 lines of shell this bundle is past the point where a self-contained template stays readable. + +| Script | Runs as | +|---|---| +| `verify_toolchain.sh` | `BioToolchain` environment `onEnter` | +| `build_index.sh` | `BuildIndex` | +| `qc_sample.sh` | `QualityControl` | +| `align_sample.sh` | `AlignReads` | +| `call_region.sh` | `CallVariants` | +| `merge_variants.sh` | `MergeVariants` | +| `common.sh` | sourced by the others, never run directly | + +A hidden `JobScriptDir` PATH parameter stages the directory. The tradeoff is that scripts in it +cannot use `{{Param.Name}}` substitution, since they are ordinary files rather than templated +embedded files. Values arrive as `--flag=value` arguments instead: + +```yaml +args: +- '{{Param.JobScriptDir}}/call_region.sh' +- '--regions={{Param.Regions}}' +- '--region-index={{Task.Param.RegionIndex}}' +``` + +Each value is labeled rather than positional, because `call_region.sh` takes ten of them and an +unlabeled list of ten is easy to reorder by accident. Each script rejects an unrecognized flag and +names any required flag it did not receive, so a typo fails on the first task instead of silently +taking a default. + +`common.sh` exists for one reason worth knowing: `call_region.sh` writes the per-region VCFs and +`merge_variants.sh` looks for them, so the two must derive the same filename. That derivation is +`region_vcf_path`, defined once and sourced by both, rather than duplicated with a comment asking +each copy to stay in step. `parse_list` is there for the same reason, since every script splits a +comma-separated parameter the same way. + +Editing a script does not require touching the template, which also means `openjd run` picks up a +change with no re-submission. + +### How steps pass files to each other + +Each step runs in its own session on a worker, with its own session directory. That shapes how this +bundle is written, in two ways. + +**Every path is derived from a job parameter, never passed between steps.** A path one step writes +down does not exist for the next, whose session directory is elsewhere. `BuildIndex` collects the +reference, its `.fai`, and the bwa index into `OutputDir/reference/`, and the later steps recompute +that location from the `OutputDir` and `ReferenceFastaName` parameters, which each session resolves +for itself. + +**A `PATH` parameter stages exactly what it names.** Job attachments uploads the single path a +`FILE`-typed parameter points at and does not sweep in files beside it. That is why the reference is a +`ReferenceDir` directory plus a `ReferenceFastaName` filename rather than one `ReferenceFasta` file +parameter: a `.fai` or bwa index sitting next to the FASTA has to be inside a staged *directory* to +reach the worker at all, and without that, `BuildIndex`'s reuse branches could never fire on a real +submission. + +**A step must declare every step whose output it reads.** A step's `dependencies` entitle it to those +steps' outputs, and the entitlement does not chain: given `A → B → C`, step `C` is not promised `A`'s +outputs without its own `dependsOn: A`. `CallVariants` lists both `AlignReads` and `BuildIndex`, +since it needs the BAMs as well as the reference index. `MergeVariants` lists all four earlier steps, +because its MultiQC report covers FastQC output and flagstat summaries as well as the VCFs it merges. + +That entitlement is a lower bound: job attachments guarantees the declared inputs and depended-on +outputs are present, not that nothing else is. Extra files may well be there, since a reused session +leaves earlier state behind and one edge brings everything that step wrote. A step reading an +undeclared file can succeed by luck and fail later when scheduling differs, so declare what you read. + +### Why the `...Range` parameters exist + +A `jobtemplate-2023-09` task parameter range cannot be computed from the length of another parameter, +so `Samples` and `Regions` are each paired with a range that sizes the parameter space and has to be +kept in step with its list. + +Each range must cover its whole list, so `SampleRange` and `RegionRange` are `0-1` and `0-3` to match +the default lists. Change a list and its range together. + +A range reaching **past** the end of its list is safe: the extra tasks detect it and exit without +doing work. A range that **skips** an entry is not, and fails rather than producing a partial result: +`CallVariants` requires a BAM for every entry in `Samples` and `MergeVariants` a VCF for every entry +in `Regions`, and both name the missing entries and the range that would cover them. For the same +reason neither step globs the output directory for its inputs. Both build the expected filenames, so +files left by an earlier run cannot quietly join the result. + +**To work on a subset, shorten the list rather than the range.** `OutputDir` is declared `dataFlow: +OUT`, so job attachments treats it as output only and a new job starts with it empty. A narrowed range +has no earlier results to combine with, so the gather insists on a complete set. Calling one region +means `-p Regions=1:137000-139999 -p RegionRange=0`, not `-p RegionRange=2`: + +```console +deadline bundle submit . \ + -p Samples=tiny_n -p SampleRange=0 \ + -p Regions=1:137000-139999 -p RegionRange=0 +``` + +Supporting a true incremental rerun (keeping earlier per-region VCFs and adding to them) needs +`OutputDir` declared `INOUT` so the prior results are staged back in. That is a deliberate design +choice rather than an oversight: `INOUT` uploads and re-downloads the whole output directory on every +run, and for a sample the simpler contract is worth more than incremental reruns. + +A task parameter range is capped at 1024 elements, which is the real ceiling on the fan-out of either +axis. + +### Scaling the region list past the parameter limit + +A job parameter string is capped at 1024 characters, so `Regions` holds roughly 100 windows before it +runs out of room. A real whole-genome scatter wants more: 3-Mb windows over GRCh38 come to about a +thousand, which is also where the 1024-element range cap lands. + +To go past that, keep the list in a file rather than a parameter. Add a `FILE`-typed `PATH` parameter +beside the existing ones, one region per line: + +```yaml +- name: RegionsFile + type: PATH + objectType: FILE + dataFlow: IN + default: sample_inputs/regions.txt + description: > + One region per line, used instead of the Regions parameter when set. Lifts the + roughly 100-entry ceiling that the 1024-character parameter limit imposes. + userInterface: + control: CHOOSE_INPUT_FILE + label: Regions File + groupLabel: Parallelism +``` + +Pass it to `call_region.sh` and `merge_variants.sh` as `--regions-file={{Param.RegionsFile}}`, and read +it in `common.sh` alongside the existing `parse_list`, so both steps still derive region names and the +per-region VCF filenames from one implementation: + +```bash +# Read one region per line, ignoring blank lines and '#' comments. +read_region_file() { + local -n _out="$1" + local _path="$2" _line + _out=() + while IFS= read -r _line || [[ -n "$_line" ]]; do + _line="${_line%%#*}" + _line="${_line#"${_line%%[![:space:]]*}"}" + _line="${_line%"${_line##*[![:space:]]}"}" + [[ -n "$_line" ]] && _out+=("$_line") + done < "$_path" +} +``` + +`RegionRange` still has to cover the file's line count, so the range stays a parameter that must be +kept in step. Generating both from a `.fai` is the usual approach, since `cut -f1,2 .fai` +gives the contig lengths that windowing needs: + +```console +# Write 3-Mb windows and report the range that covers them +awk 'BEGIN{OFS=""} {for (s=1; s<=$2; s+=3000000) {e=s+2999999; if (e>$2) e=$2; print $1,":",s,"-",e}}' \ + reference/human_g1k_v37_decoy.small.fasta.fai > sample_inputs/regions.txt +echo "RegionRange=0-$(( $(wc -l < sample_inputs/regions.txt) - 1 ))" +``` + +This bundle keeps the parameter form because a reader can see the whole scatter in the submission +command, and the sample data needs four windows. The file form is the change to make when the region +count outgrows what a parameter can carry. + +The EXPR extension's +[`LIST[STRING]`](https://github.com/OpenJobDescription/openjd-specifications/wiki/2023-09-Template-Schemas#211-jobliststringparameterdefinition-extension-expr) +type removes the need for the paired range, at the cost of requiring an implementation that supports +the extension. No other sample here uses EXPR yet, so this bundle stays on the base 2023-09 schema. + +## Security, cost, and cleanup + +Running this job incurs Deadline Cloud worker, storage, and data transfer charges. + +`fetch_test_data.py` downloads from a public, unauthenticated S3 bucket over HTTPS and writes only +inside `sample_inputs/`. Re-running it skips files that already exist. `--clean` removes them first. + +Genomic sequence data is often subject to consent, privacy, and jurisdictional restrictions. Set +queue, S3, and farm access controls to match before pointing this bundle at real data, and prefer a +dedicated queue whose job attachments bucket has the encryption and access logging your data +governance requires. + +The `output/` directory is excluded by this repository's `.gitignore`. Delete it to reclaim space. + +## Troubleshooting + +**The job fails immediately with "required tool ... is not on PATH."** The queue has no conda +queue environment, or its channels or packages are wrong. The error lists the exact values to set. + +**`bcftools mpileup` reports an unknown region, or a region yields zero variants.** Region names +must match the reference exactly. GRCh37-style references name contigs `1`, `2`, `X`. GRCh38 +references from UCSC and the GATK resource bundle name them `chr1`, `chr2`, `chrX`. Check with +`samtools faidx && cut -f1 .fai`. + +**"no VCF for N of M requested region(s)" or "no alignment for N of M requested sample(s)."** +`RegionRange` or `SampleRange` does not cover its whole list, so some entries were never processed. +Both errors name the missing entries and the range that would cover them. If you narrowed a range +deliberately to redo part of a run, shorten the matching list instead. See +[Why the `...Range` parameters exist](#why-the-range-parameters-exist). + +**A sample is missing from the VCF's sample columns.** `AlignReads` writes the sample name into the +BAM's `@RG` line, which is what `bcftools` reads. If you supply pre-aligned BAMs from another +pipeline, confirm they carry an `@RG` line with `SM:` set. + +**Fewer variants than expected.** `MinQual` (default 20) and `MinDepth` (default 5) discard +low-confidence calls. On shallow data those defaults can remove most sites; set both to 0 to see the +unfiltered output, and read the counts `CallVariants` logs per region. + +**`bcftools index` fails with a compression error.** The index requires BGZF-compressed input, not +plain gzip. `bcftools call --output-type z` produces BGZF, so this only arises if you substitute +your own compression step. + +**A reference contig longer than 512 Mb fails to index.** Both `CallVariants` and `MergeVariants` +index with `--tbi`, whose format caps contig length at 2^29 bases. Every human chromosome fits +comfortably, but some plant and amphibian genomes do not; switch those to CSI by replacing `--tbi` +with `--csi` in both steps. + +Worker logs for a failed task are available through the Deadline Cloud monitor, or with +`deadline job logs --job-id `, which also accepts `--session-id` to narrow the output to one +step. See the +[log retrieval guidance](https://docs.aws.amazon.com/deadline-cloud/latest/userguide/view-logs.html). + +While iterating, submit with `--max-retries-per-task 0` so a broken step fails once instead of +retrying before the job gives up. + +## Related job bundles + +* [ESMFold protein structure prediction](../esmfold_predict/): GPU bioinformatics, FASTA to PDB +* [GROMACS molecular dynamics](../gromacs_md/): multi-stage simulation with replica fan-out +* [AutoDock Vina virtual screening](../virtual_screening_vina/): chunked molecular docking +* [Monte Carlo simulation](../monte_carlo_simulation/): the fan-out/fan-in pattern with task chunking + +## Related resources + +* [bioconda](https://bioconda.github.io/): channel setup and the package index +* [Open Job Description step parameter space definitions](https://github.com/OpenJobDescription/openjd-specifications/wiki/2023-09-Template-Schemas#34-stepparameterspacedefinition) diff --git a/job_bundles/variant_calling_bwa/sample_inputs/.gitignore b/job_bundles/variant_calling_bwa/sample_inputs/.gitignore new file mode 100644 index 00000000..8bb70cef --- /dev/null +++ b/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/ diff --git a/job_bundles/variant_calling_bwa/sample_inputs/README.md b/job_bundles/variant_calling_bwa/sample_inputs/README.md new file mode 100644 index 00000000..d4c14d13 --- /dev/null +++ b/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. diff --git a/job_bundles/variant_calling_bwa/sample_inputs/fetch_test_data.py b/job_bundles/variant_calling_bwa/sample_inputs/fetch_test_data.py new file mode 100644 index 00000000..bbf14895 --- /dev/null +++ b/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 ``_R1.fastq.gz`` / +``_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()) diff --git a/job_bundles/variant_calling_bwa/scripts/align_sample.sh b/job_bundles/variant_calling_bwa/scripts/align_sample.sh new file mode 100644 index 00000000..f1e50e91 --- /dev/null +++ b/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 +# ".bwt" or the ".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" diff --git a/job_bundles/variant_calling_bwa/scripts/build_index.sh b/job_bundles/variant_calling_bwa/scripts/build_index.sh new file mode 100644 index 00000000..590c9c14 --- /dev/null +++ b/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 ".64.*" instead of ".*", 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" diff --git a/job_bundles/variant_calling_bwa/scripts/call_region.sh b/job_bundles/variant_calling_bwa/scripts/call_region.sh new file mode 100644 index 00000000..dd05eced --- /dev/null +++ b/job_bundles/variant_calling_bwa/scripts/call_region.sh @@ -0,0 +1,105 @@ +#!/usr/bin/env bash +# Call variants over one reference region, jointly across all samples. +set -euo pipefail +source "$(dirname "${BASH_SOURCE[0]}")/common.sh" + +REGIONS_CSV="" +REGION_INDEX="" +SAMPLES_CSV="" +OUT="" +REFERENCE_NAME="" +MIN_MAPPING_QUALITY="" +MIN_BASE_QUALITY="" +PLOIDY="" +MIN_QUAL="" +MIN_DEPTH="" +for arg in "$@"; do + case "$arg" in + --regions=*) REGIONS_CSV="${arg#*=}" ;; + --region-index=*) REGION_INDEX="${arg#*=}" ;; + --samples=*) SAMPLES_CSV="${arg#*=}" ;; + --output-dir=*) OUT="${arg#*=}" ;; + --reference-name=*) REFERENCE_NAME="${arg#*=}" ;; + --min-mapping-quality=*) MIN_MAPPING_QUALITY="${arg#*=}" ;; + --min-base-quality=*) MIN_BASE_QUALITY="${arg#*=}" ;; + --ploidy=*) PLOIDY="${arg#*=}" ;; + --min-qual=*) MIN_QUAL="${arg#*=}" ;; + --min-depth=*) MIN_DEPTH="${arg#*=}" ;; + *) unknown_arg "$arg" ;; + esac +done +require_arg "$REGIONS_CSV" --regions +require_arg "$REGION_INDEX" --region-index +require_arg "$SAMPLES_CSV" --samples +require_arg "$OUT" --output-dir +require_arg "$REFERENCE_NAME" --reference-name +require_arg "$MIN_MAPPING_QUALITY" --min-mapping-quality +require_arg "$MIN_BASE_QUALITY" --min-base-quality +require_arg "$PLOIDY" --ploidy +require_arg "$MIN_QUAL" --min-qual +require_arg "$MIN_DEPTH" --min-depth + +parse_list REGIONS "$REGIONS_CSV" +require_non_empty_list "${#REGIONS[@]}" Regions +parse_list SAMPLES "$SAMPLES_CSV" +require_non_empty_list "${#SAMPLES[@]}" Samples +validate_sample_names "${SAMPLES[@]}" +if (( REGION_INDEX >= ${#REGIONS[@]} )); then + echo "Region index $REGION_INDEX is past the end of the region list; nothing to do." + exit 0 +fi +REGION="${REGIONS[$REGION_INDEX]}" +BAM_DIR="$OUT/alignments" +VCF_DIR="$OUT/vcf_by_region" +mkdir -p "$VCF_DIR" + +VCF="$(region_vcf_path "$VCF_DIR" "$REGION_INDEX" "$REGION")" + +# mpileup requires a .fai beside the reference it is given, so use the copy +# BuildIndex assembled. Derived from the output directory, not read from a file. +REF="$OUT/reference/$REFERENCE_NAME" +if [[ ! -f "${REF}.fai" ]]; then + echo "ERROR: no .fai index at ${REF}.fai. Did BuildIndex run?" >&2 + exit 1 +fi + +# Build the BAM list from the sample list rather than globbing the directory. A +# glob would silently pick up BAMs left over from an earlier run or a reused +# session and call a cohort the user did not ask for. +BAMS=() +for s in "${SAMPLES[@]}"; do + bam="$BAM_DIR/${s}.sorted.bam" + if [[ ! -f "$bam" ]]; then + echo "ERROR: no alignment for sample '$s' at $bam." >&2 + echo "AlignReads must run for every sample listed in Samples." >&2 + exit 1 + fi + BAMS+=("$bam") +done +echo "Calling region '$REGION' jointly across ${#BAMS[@]} sample(s): ${SAMPLES[*]}" + +# --ploidy takes an assembly preset; '1' and '2' are the presets meaning treat +# every sample as haploid or diploid respectively. +bcftools mpileup \ + --fasta-ref "$REF" \ + --regions "$REGION" \ + -q "$MIN_MAPPING_QUALITY" \ + -Q "$MIN_BASE_QUALITY" \ + --annotate FORMAT/AD,FORMAT/DP \ + --output-type u \ + "${BAMS[@]}" \ + | bcftools call \ + --multiallelic-caller \ + --variants-only \ + --ploidy "$PLOIDY" \ + --output-type u \ + | bcftools filter \ + --include "QUAL>=$MIN_QUAL && INFO/DP>=$MIN_DEPTH" \ + --output-type z \ + --output "$VCF" + +# An index is required here, not just convenient: MergeVariants uses +# 'bcftools concat --allow-overlaps', which needs indexed inputs. +bcftools index --force --tbi "$VCF" +COUNT="$(bcftools view --no-header "$VCF" | wc -l)" +echo "Region '$REGION': ${COUNT} variant(s) passing QUAL>=$MIN_QUAL DP>=$MIN_DEPTH -> $VCF" diff --git a/job_bundles/variant_calling_bwa/scripts/common.sh b/job_bundles/variant_calling_bwa/scripts/common.sh new file mode 100644 index 00000000..d7fbf0ac --- /dev/null +++ b/job_bundles/variant_calling_bwa/scripts/common.sh @@ -0,0 +1,108 @@ +#!/usr/bin/env bash +# Helpers shared by the bundled scripts. Source it, do not run it: +# +# source "$(dirname "${BASH_SOURCE[0]}")/common.sh" +# +# Nothing here calls the bioinformatics tools. This file exists so that the two +# scripts which have to agree on a filename derive it from one implementation +# instead of two comments asking each other to stay in step. + +# Split a comma-separated job parameter into the named array, trimming +# surrounding whitespace from every element so that "a, b" does not yield an +# entry with a leading space, which would otherwise surface as a confusing +# missing-file error or a filename containing a space. +# +# Only the surrounding whitespace goes. Deleting interior spaces too would +# rewrite 'tiny n' to 'tinyn' before validate_sample_names ever sees it, and the +# job would then fail on a filename the user never typed rather than on the name +# they did. +# +# parse_list SAMPLES "tiny_n, tiny_t" +parse_list() { + local -n _out="$1" + local _raw="$2" + local _i _elem + IFS=',' read -r -a _out <<< "$_raw" + for _i in "${!_out[@]}"; do + _elem="${_out[$_i]}" + _elem="${_elem#"${_elem%%[![:space:]]*}"}" + _elem="${_elem%"${_elem##*[![:space:]]}"}" + _out[$_i]="$_elem" + done +} + +# Reject a list parameter that has no entries. +# +# An empty parameter parses to zero elements, which every per-index task reports +# as "past the end of the list" while the gather's completeness checks find +# nothing missing to complain about. Without this the job would fail deep inside +# a tool -- bcftools handed no BAM arguments -- instead of naming the parameter +# that was left empty. +# +# require_non_empty_list "${#SAMPLES[@]}" Samples +require_non_empty_list() { + local _count="$1" _param="$2" + if (( _count == 0 )); then + echo "ERROR: the $_param parameter is empty; it needs at least one entry." >&2 + exit 1 + fi +} + +# Sample names are pasted into input, output, and temporary paths, so reject +# anything that is not a plain filename component. A name like 'a/../../b' would +# otherwise resolve outside the directory it looks like it is in. +# +# Letters, digits, dot, underscore, and hyphen cover the sample naming that +# sequencing platforms produce; a leading dot is refused so a name cannot be +# '..' or produce a hidden file. +validate_sample_names() { + local _name + for _name in "$@"; do + if [[ -z "$_name" ]]; then + echo "ERROR: empty sample name in the Samples list." >&2 + exit 1 + fi + if [[ ! "$_name" =~ ^[A-Za-z0-9_][A-Za-z0-9._-]*$ ]]; then + echo "ERROR: invalid sample name '$_name'." >&2 + echo "Sample names may contain letters, digits, '.', '_', and '-', and must" >&2 + echo "not start with '.', because they are used to build file paths." >&2 + exit 1 + fi + done +} + +# Region strings may contain ':' and '*', which are not safe in filenames. +safe_region_name() { + printf '%s' "$1" | tr ':*/|' '____' +} + +# The per-region VCF path. CallVariants writes it and MergeVariants looks for +# it, so both must agree exactly; that is why this lives here rather than being +# spelled out twice. +# +# The index prefix keeps the name unique: the character mapping above is not +# injective ('1:100-200' and '1|100-200' both become '1_100-200'), and two +# parallel CallVariants tasks writing one file would corrupt it silently. +# +# region_vcf_path "$OUT/vcf_by_region" 2 "1:134000-136999" +region_vcf_path() { + local _dir="$1" _index="$2" _region="$3" + printf '%s/region_%04d_%s.vcf.gz' \ + "$_dir" "$_index" "$(safe_region_name "$_region")" +} + +# Fail with a message naming the flag a script needs but did not receive. +require_arg() { + local _value="$1" _flag="$2" + if [[ -z "$_value" ]]; then + echo "ERROR: missing required argument $_flag" >&2 + exit 1 + fi +} + +# Reject an unrecognized flag rather than ignoring it, so a typo in the job +# template fails loudly on the first task instead of silently using a default. +unknown_arg() { + echo "ERROR: unrecognized argument: $1" >&2 + exit 1 +} diff --git a/job_bundles/variant_calling_bwa/scripts/merge_variants.sh b/job_bundles/variant_calling_bwa/scripts/merge_variants.sh new file mode 100644 index 00000000..b47df340 --- /dev/null +++ b/job_bundles/variant_calling_bwa/scripts/merge_variants.sh @@ -0,0 +1,114 @@ +#!/usr/bin/env bash +# Gather the per-region VCFs into one sorted, indexed VCF, summarize it, and +# build a MultiQC report over the run. +set -euo pipefail +source "$(dirname "${BASH_SOURCE[0]}")/common.sh" + +REGIONS_CSV="" +SAMPLES_CSV="" +OUT="" +REFERENCE_NAME="" +for arg in "$@"; do + case "$arg" in + --regions=*) REGIONS_CSV="${arg#*=}" ;; + --samples=*) SAMPLES_CSV="${arg#*=}" ;; + --output-dir=*) OUT="${arg#*=}" ;; + --reference-name=*) REFERENCE_NAME="${arg#*=}" ;; + *) unknown_arg "$arg" ;; + esac +done +require_arg "$REGIONS_CSV" --regions +require_arg "$SAMPLES_CSV" --samples +require_arg "$OUT" --output-dir +require_arg "$REFERENCE_NAME" --reference-name + +VCF_DIR="$OUT/vcf_by_region" +MERGED="$OUT/variants.vcf.gz" + +# Collect exactly the VCFs the requested regions should have produced, in the +# order the regions were given. Globbing the directory instead would pick up +# stale files from an earlier run with different regions and merge a result the +# user did not ask for. +parse_list REGIONS "$REGIONS_CSV" +require_non_empty_list "${#REGIONS[@]}" Regions +VCFS=() +MISSING=() +for i in "${!REGIONS[@]}"; do + vcf="$(region_vcf_path "$VCF_DIR" "$i" "${REGIONS[$i]}")" + if [[ -f "$vcf" ]]; then + VCFS+=("$vcf") + else + MISSING+=("${REGIONS[$i]}") + fi +done + +# A missing VCF means CallVariants never ran for that region, which happens when +# RegionRange does not cover the whole region list. +if (( ${#MISSING[@]} > 0 )); then + echo "ERROR: no VCF for ${#MISSING[@]} of ${#REGIONS[@]} requested region(s):" >&2 + printf ' %s\n' "${MISSING[@]}" >&2 + echo "RegionRange must cover every region: '0-$(( ${#REGIONS[@]} - 1 ))' for this list." >&2 + exit 1 +fi +# Same check for samples: a SampleRange that skips an entry would leave that +# sample unaligned and silently absent from the VCF. +parse_list SAMPLES "$SAMPLES_CSV" +require_non_empty_list "${#SAMPLES[@]}" Samples +validate_sample_names "${SAMPLES[@]}" +MISSING_BAM=() +for s in "${SAMPLES[@]}"; do + [[ -f "$OUT/alignments/${s}.sorted.bam" ]] || MISSING_BAM+=("$s") +done +if (( ${#MISSING_BAM[@]} > 0 )); then + echo "ERROR: no alignment for ${#MISSING_BAM[@]} of ${#SAMPLES[@]} requested sample(s):" >&2 + printf ' %s\n' "${MISSING_BAM[@]}" >&2 + echo "SampleRange must cover every sample: '0-$(( ${#SAMPLES[@]} - 1 ))' for this list." >&2 + exit 1 +fi +echo "Concatenating ${#VCFS[@]} per-region VCF(s)." + +# --allow-overlaps lets concat accept inputs that are not already in coordinate +# order, which matters because the region list can be in any order. --rm-dups +# exact then drops records that appear in more than one input: -a alone tolerates +# overlap but does not deduplicate, so padded or overlapping regions would +# otherwise produce repeated records. +# +# bcftools norm left-aligns indels and splits multiallelic records, which is what +# makes the VCF comparable against another callset. It runs after concat because +# left-alignment can move an indel's position. +REF="$OUT/reference/$REFERENCE_NAME" +bcftools concat --allow-overlaps --rm-dups exact --output-type u "${VCFS[@]}" \ + | bcftools norm --fasta-ref "$REF" --multiallelics -any \ + --output-type z --output "$MERGED" +bcftools index --force --tbi "$MERGED" + +SUMMARY="$OUT/variant_summary.txt" +{ + echo "=== Variant Calling Summary ===" + echo "Regions called: ${#VCFS[@]}" + echo "Total variants: $(bcftools view --no-header "$MERGED" | wc -l)" + echo "Samples: $(bcftools view -h "$MERGED" | grep '^#CHROM' | cut -f10- | tr '\t' ' ')" + echo "" + echo "--- variants per contig ---" + bcftools view --no-header "$MERGED" | awk '{print $1}' | sort | uniq -c \ + | awk '{printf " %-12s %s\n", $2, $1}' + echo "" + echo "--- bcftools stats ---" + # '|| true' on the greps: a no-match grep is exit 1, and as the last command of + # this redirected group that would abort the script leaving $SUMMARY truncated + # with nothing printed, since stdout is captured. + bcftools stats "$MERGED" | { grep -E '^SN' || true; } | cut -f3- +} > "$SUMMARY" +cat "$SUMMARY" + +# MultiQC exits 1 when it finds no analysis results to summarize, which is not a +# pipeline failure: the VCF above is the actual result, and the report is a +# convenience over it. Keep its output visible so a genuine error stays +# diagnosable. +echo "Building MultiQC report..." +if multiqc --force --outdir "$OUT/multiqc" "$OUT"; then + echo "MultiQC report: $OUT/multiqc/multiqc_report.html" +else + echo "MultiQC found nothing to summarize (exit $?); the merged VCF is unaffected." +fi +echo "Merged VCF: $MERGED" diff --git a/job_bundles/variant_calling_bwa/scripts/qc_sample.sh b/job_bundles/variant_calling_bwa/scripts/qc_sample.sh new file mode 100644 index 00000000..0707df07 --- /dev/null +++ b/job_bundles/variant_calling_bwa/scripts/qc_sample.sh @@ -0,0 +1,54 @@ +#!/usr/bin/env bash +# Run FastQC on one sample's read pair. One task per sample. +set -euo pipefail +source "$(dirname "${BASH_SOURCE[0]}")/common.sh" + +SAMPLES_CSV="" +SAMPLE_INDEX="" +READS="" +OUT="" +for arg in "$@"; do + case "$arg" in + --samples=*) SAMPLES_CSV="${arg#*=}" ;; + --sample-index=*) SAMPLE_INDEX="${arg#*=}" ;; + --reads-dir=*) READS="${arg#*=}" ;; + --output-dir=*) OUT="${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 + +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]}" +QC_DIR="$OUT/qc" +mkdir -p "$QC_DIR" + +R1="$READS/${SAMPLE}_R1.fastq.gz" +R2="$READS/${SAMPLE}_R2.fastq.gz" +for f in "$R1" "$R2"; do + if [[ ! -f "$f" ]]; then + echo "ERROR: FASTQ not found: $f" >&2 + exit 1 + fi +done + +echo "Running FastQC on $SAMPLE" +fastqc --outdir "$QC_DIR" --threads 2 "$R1" "$R2" + +# Name the reports this task produced, so the log says what to look for in the +# output directory. Derived from the two input filenames rather than globbing +# $QC_DIR, which is shared: a glob would also list the reports left by the other +# samples' tasks when a session is reused. +for f in "$R1" "$R2"; do + echo " report: $QC_DIR/$(basename "$f" .fastq.gz)_fastqc.html" +done +echo "QC complete for $SAMPLE" diff --git a/job_bundles/variant_calling_bwa/scripts/verify_toolchain.sh b/job_bundles/variant_calling_bwa/scripts/verify_toolchain.sh new file mode 100644 index 00000000..3674eb19 --- /dev/null +++ b/job_bundles/variant_calling_bwa/scripts/verify_toolchain.sh @@ -0,0 +1,49 @@ +#!/usr/bin/env bash +# Verify the tools every step needs are on PATH before any of them runs. +# Runs as the job environment's onEnter action. +set -euo pipefail +source "$(dirname "${BASH_SOURCE[0]}")/common.sh" + +REQUIRED_TOOLS_CSV="" +CONDA_PACKAGES="" +CONDA_CHANNELS="" +for arg in "$@"; do + case "$arg" in + --required-tools=*) REQUIRED_TOOLS_CSV="${arg#*=}" ;; + --conda-packages=*) CONDA_PACKAGES="${arg#*=}" ;; + --conda-channels=*) CONDA_CHANNELS="${arg#*=}" ;; + *) unknown_arg "$arg" ;; + esac +done +require_arg "$REQUIRED_TOOLS_CSV" --required-tools + +# Checking here means a queue environment missing a tool fails once, up front, +# with the packages and channels to fix it, rather than partway into a step. +parse_list REQUIRED_TOOLS "$REQUIRED_TOOLS_CSV" +echo "Verifying toolchain: ${REQUIRED_TOOLS[*]}" +missing=0 +for tool in "${REQUIRED_TOOLS[@]}"; do + if ! command -v "$tool" >/dev/null 2>&1; then + echo "ERROR: required tool '$tool' is not on PATH." >&2 + missing=1 + fi +done +if [[ "$missing" -ne 0 ]]; then + echo "" >&2 + echo "Attach a conda queue environment that provides:" >&2 + echo " packages: $CONDA_PACKAGES" >&2 + echo " channels: $CONDA_CHANNELS" >&2 + exit 1 +fi + +# Log the htslib tool versions, since those are the ones whose behavior differs +# enough between releases to matter when reproducing a callset. Not every tool in +# the list, because they do not agree on a version flag. +# +# Read the first line without a pipe: 'set -o pipefail' plus 'head' closing the +# pipe early makes "tool | head" exit 141 (SIGPIPE) even though the tool +# succeeded. +for tool in samtools bcftools; do + read -r line < <("$tool" --version) && echo "$line" +done +echo "Toolchain ready." diff --git a/job_bundles/variant_calling_bwa/template.yaml b/job_bundles/variant_calling_bwa/template.yaml new file mode 100644 index 00000000..035b618c --- /dev/null +++ b/job_bundles/variant_calling_bwa/template.yaml @@ -0,0 +1,459 @@ +specificationVersion: 'jobtemplate-2023-09' +name: Genetic Variant Calling from DNA Sequencing +description: | + Finds the genetic differences between an individual and a reference genome. + + A sequencing machine produces hundreds of millions of short DNA fragments from + random positions rather than reading a genome end to end. This job aligns those + fragments to a reference to work out where each came from (bwa), compares the + stacked-up fragments against the reference to decide where the individual + genuinely differs (bcftools), and merges the results into one list of variants. + Input is FASTQ, output is VCF. "Germline" means the variants an individual was + born with, rather than ones a tumor acquired later. + + Both expensive stages divide cleanly: each sample is aligned on its own worker, + and separate regions of the genome are then called in parallel (the scatter) + and stitched back together (the gather). + + The tool sequence follows the AWS HealthOmics WDL variant-calling tutorial + pipeline, reimplemented for Open Job Description: + https://github.com/aws-samples/aws-healthomics-tutorials/tree/main/example-workflows/wdl/variant-calling-pipeline + + Two things differ. That pipeline scatters over samples and then gathers every + sample into one whole-genome call; this bundle adds a second scatter over + regions so calling itself is distributed. It also takes paired-end reads, which + sequencers usually produce in pairs, where the WDL takes one FASTQ per sample. + + Where that pipeline runs each tool in its own container image, this bundle + declares them as conda packages. Those images are BioContainers, which are + themselves built from bioconda, so the tool builds are equivalent. + + Requires a Deadline Cloud queue with a conda queue environment whose channels + include conda-forge and bioconda. See ../../queue_environments/. + +parameterDefinitions: + +# --- Input / Output --- +- name: ReadsDir + type: PATH + objectType: DIRECTORY + dataFlow: IN + default: sample_inputs/reads + description: > + Directory holding paired-end FASTQ files. For each sample name S, this + directory must contain S_R1.fastq.gz and S_R2.fastq.gz. + userInterface: + control: CHOOSE_DIRECTORY + label: Reads Directory + groupLabel: Input / Output + +- name: ReferenceDir + type: PATH + objectType: DIRECTORY + dataFlow: IN + default: sample_inputs/reference + description: > + Directory holding the reference genome FASTA and any indexes already built + beside it. BuildIndex reuses a staged .fai and bwa index when it finds them + and builds whichever are missing. + userInterface: + control: CHOOSE_DIRECTORY + label: Reference Directory + groupLabel: Input / Output + +- name: ReferenceFastaName + type: STRING + default: human_g1k_v37_decoy.small.fasta + description: > + Filename of the reference FASTA within the reference directory. A name only, + not a path, so that a directory holding several genomes selects one. + userInterface: + control: LINE_EDIT + label: Reference FASTA Name + groupLabel: Input / Output + +- name: OutputDir + type: PATH + objectType: DIRECTORY + dataFlow: OUT + default: output + description: "Directory for alignments, per-region VCFs, the merged VCF, and QC reports." + userInterface: + control: CHOOSE_DIRECTORY + label: Output Directory + groupLabel: Input / Output + +# --- Parallelism --- +# +# Samples and Regions are comma-separated lists, each paired with a range that +# selects which entries to process. A 2023-09 task parameter range cannot be +# computed from another parameter's length, so the range is a parameter of its +# own and must be kept in step with its list. The scripts index into the list +# and no-op when the index is past the end, so a range reaching beyond the list +# is safe; MergeVariants fails loudly if any listed entry was skipped, since a +# range that omits one would otherwise drop it from the result silently. +# The EXPR extension's LIST[STRING] type removes the need for the paired range. +- name: Samples + type: STRING + default: "tiny_n,tiny_t" + description: > + Comma-separated sample names to align in parallel. Each name S requires + S_R1.fastq.gz and S_R2.fastq.gz in the reads directory. + userInterface: + control: LINE_EDIT + label: Samples + groupLabel: Parallelism + +- name: SampleRange + type: STRING + default: "0-1" + description: > + Which entries of Samples to align, as indices from 0, so '0-1' for two + samples. This must cover the whole list: MergeVariants requires an alignment + for every name in Samples, and a job starts with an empty output directory, + so a run that skips an entry has nothing to fall back on. To work on a subset, + shorten Samples itself rather than narrowing this range. + userInterface: + control: LINE_EDIT + label: Sample Range + groupLabel: Parallelism + +- name: Regions + type: STRING + default: "1:131000-133999,1:134000-136999,1:137000-139999,1:140000-142999" + description: > + Comma-separated reference regions to call variants over, in parallel. Each + value is passed to 'bcftools mpileup -r', so it may be a whole contig (1) or + a subrange (1:131000-133999). Names must match the reference FASTA. + + The default splits contig 1 into four windows over the stretch the sample + reads align to, so the scatter tasks have real data to work on. Windowing is + also how you scatter a real genome: chromosomes vary widely in size, so one + task would run far longer than the rest. Keep the windows non-overlapping, + since overlapping ones would call the same site twice. + + Whichever you choose, RegionRange must cover every entry. The practical + ceiling is around 100 entries, because a job parameter string is limited to + 1024 characters. Past that, stage the list as a FILE-typed PATH parameter + holding one region per line; the README section "Scaling the region list past + the parameter limit" sketches the change. + userInterface: + control: LINE_EDIT + label: Regions + groupLabel: Parallelism + +- name: RegionRange + type: STRING + default: "0-3" + description: > + Which entries of Regions to call, as indices from 0, so '0-3' for four + regions. This must cover the whole list, for the same reason as SampleRange: + MergeVariants requires a VCF for every entry in Regions. To call a subset, + shorten Regions itself rather than narrowing this range. + userInterface: + control: LINE_EDIT + label: Region Range + groupLabel: Parallelism + +# --- Calling Parameters --- +- name: MinMappingQuality + type: INT + default: 20 + minValue: 0 + maxValue: 60 + description: "Skip alignments with mapping quality below this value (bcftools mpileup -q)." + userInterface: + control: SPIN_BOX + label: Min Mapping Quality + groupLabel: Calling Parameters + +- name: MinBaseQuality + type: INT + default: 20 + minValue: 0 + maxValue: 60 + description: "Skip bases with base quality below this value (bcftools mpileup -Q)." + userInterface: + control: SPIN_BOX + label: Min Base Quality + groupLabel: Calling Parameters + +- name: Ploidy + type: STRING + default: "2" + allowedValues: ["1", "2"] + description: > + Ploidy preset for bcftools call: 2 treats every sample as diploid, 1 as + haploid. This applies to the whole job; a real analysis needs per-contig + ploidy so that chrX, chrY, and the mitochondrion are handled correctly. + userInterface: + control: DROPDOWN_LIST + label: Ploidy + groupLabel: Calling Parameters + +- name: MinQual + type: FLOAT + default: 20.0 + minValue: 0.0 + description: > + Discard called sites below this QUAL. 'bcftools call --variants-only' emits + every variant site, not only confident ones, so without a filter roughly half + the raw calls on a small test dataset are noise. Set to 0 to keep every site. + userInterface: + control: SPIN_BOX + label: Min Variant QUAL + groupLabel: Calling Parameters + +- name: MinDepth + type: INT + default: 5 + minValue: 0 + description: > + Discard called sites whose total depth (INFO/DP) is below this. Guards against + genotypes asserted from one or two reads. Set to 0 to keep every site. + userInterface: + control: SPIN_BOX + label: Min Depth + groupLabel: Calling Parameters + +# --- Software Environment --- +- name: CondaPackages + type: STRING + default: "bwa samtools bcftools fastqc multiqc" + description: > + Packages the conda queue environment installs. bwa, samtools, bcftools and + fastqc come from bioconda; multiqc is noarch. + userInterface: + control: LINE_EDIT + label: Conda Packages + groupLabel: Software Environment + +- name: CondaChannels + type: STRING + default: "conda-forge bioconda" + description: > + Conda channels, highest priority first. bioconda depends on conda-forge for + its runtime libraries and documents this order, so conda-forge must come + first for packages to resolve correctly. + userInterface: + control: LINE_EDIT + label: Conda Channels + groupLabel: Software Environment + +- name: RequiredTools + type: STRING + default: "bwa,samtools,bcftools,fastqc,multiqc" + description: > + Comma-separated commands the BioToolchain environment requires on PATH before + any step runs. Every step's tools are listed, so a queue environment missing + one fails immediately with the packages and channels to add rather than + partway through a step. Keep this in step with CondaPackages. + userInterface: + control: LINE_EDIT + label: Required Tools + groupLabel: Software Environment + +# The scripts live in a directory carried with the bundle rather than embedded in +# this template. A PATH parameter stages them, which means the scripts cannot use +# "{{Param.Name}}" substitution themselves; values reach them as named command +# line arguments instead. See scripts/common.sh for the helpers they share. +- name: JobScriptDir + type: PATH + objectType: DIRECTORY + dataFlow: IN + default: scripts + description: "Directory containing the bundled scripts this job runs." + userInterface: + control: HIDDEN + +jobEnvironments: +- name: BioToolchain + description: "Verify the bioconda-provided tools are on PATH before any step runs." + script: + actions: + onEnter: + command: bash + args: + - '{{Param.JobScriptDir}}/verify_toolchain.sh' + - '--required-tools={{Param.RequiredTools}}' + - '--conda-packages={{Param.CondaPackages}}' + - '--conda-channels={{Param.CondaChannels}}' + timeout: 60 + +steps: + +- name: BuildIndex + description: "Create the .fai and bwa indexes once, if they are not already present." + hostRequirements: + amounts: + # bwa index and samtools faidx are both single-threaded, so asking for more + # cores would reserve capacity that sits idle for the length of the build. + - name: "amount.worker.vcpu" + min: 2 + - name: "amount.worker.memory" + min: 8192 + attributes: + - name: "attr.worker.os.family" + anyOf: ["linux"] + script: + actions: + onRun: + command: bash + args: + - '{{Param.JobScriptDir}}/build_index.sh' + - '--reference-dir={{Param.ReferenceDir}}' + - '--reference-name={{Param.ReferenceFastaName}}' + - '--output-dir={{Param.OutputDir}}' + timeout: 43200 + +- name: QualityControl + description: "Run FastQC on each sample's read pair. One task per sample." + # No dependencies: this step reads only ReadsDir, a job input, so it can start + # immediately rather than waiting behind the reference index build. + parameterSpace: + taskParameterDefinitions: + - name: SampleIndex + type: INT + range: "{{Param.SampleRange}}" + hostRequirements: + amounts: + - name: "amount.worker.vcpu" + min: 2 + - name: "amount.worker.memory" + min: 4096 + attributes: + - name: "attr.worker.os.family" + anyOf: ["linux"] + script: + actions: + onRun: + command: bash + args: + - '{{Param.JobScriptDir}}/qc_sample.sh' + - '--samples={{Param.Samples}}' + - '--sample-index={{Task.Param.SampleIndex}}' + - '--reads-dir={{Param.ReadsDir}}' + - '--output-dir={{Param.OutputDir}}' + timeout: 3600 + +- name: AlignReads + description: "Align each sample with bwa mem, then sort and index the BAM. One task per sample." + dependencies: + - dependsOn: BuildIndex + parameterSpace: + taskParameterDefinitions: + - name: SampleIndex + type: INT + range: "{{Param.SampleRange}}" + hostRequirements: + amounts: + # 2 vCPU, matching the other steps, so that any worker a fleet brings up can + # run every step of this job. + # + # Alignment is the step that most benefits from more cores, so raise this and + # the thread counts in align_sample.sh together when running real cohorts on a + # fleet whose minimum is higher. + - name: "amount.worker.vcpu" + min: 2 + - name: "amount.worker.memory" + min: 16384 + # No amount.worker.disk.scratch requirement: service-managed fleets do not + # advertise that capability, so requiring it makes every task NOT_COMPATIBLE. + # samtools sort still needs room for temp chunks -- tens of GiB for a + # whole-genome BAM -- so size the fleet's root volume accordingly. + attributes: + - name: "attr.worker.os.family" + anyOf: ["linux"] + script: + actions: + onRun: + command: bash + args: + - '{{Param.JobScriptDir}}/align_sample.sh' + - '--samples={{Param.Samples}}' + - '--sample-index={{Task.Param.SampleIndex}}' + - '--reads-dir={{Param.ReadsDir}}' + - '--output-dir={{Param.OutputDir}}' + - '--reference-name={{Param.ReferenceFastaName}}' + - '--session-work-dir={{Session.WorkingDirectory}}' + timeout: 86400 + cancelation: + mode: NOTIFY_THEN_TERMINATE + notifyPeriodInSeconds: 30 + +- name: CallVariants + description: > + Call variants over one reference region, jointly across all samples. This is + the scatter: one task per region, all able to run in parallel. + dependencies: + # Dependencies do not chain, so BuildIndex is required here even though + # AlignReads already depends on it: mpileup reads its reference and .fai. + - dependsOn: BuildIndex + - dependsOn: AlignReads + parameterSpace: + taskParameterDefinitions: + - name: RegionIndex + type: INT + range: "{{Param.RegionRange}}" + hostRequirements: + amounts: + - name: "amount.worker.vcpu" + min: 2 + - name: "amount.worker.memory" + min: 8192 + attributes: + - name: "attr.worker.os.family" + anyOf: ["linux"] + script: + actions: + onRun: + command: bash + args: + - '{{Param.JobScriptDir}}/call_region.sh' + - '--regions={{Param.Regions}}' + - '--region-index={{Task.Param.RegionIndex}}' + - '--samples={{Param.Samples}}' + - '--output-dir={{Param.OutputDir}}' + - '--reference-name={{Param.ReferenceFastaName}}' + - '--min-mapping-quality={{Param.MinMappingQuality}}' + - '--min-base-quality={{Param.MinBaseQuality}}' + - '--ploidy={{Param.Ploidy}}' + - '--min-qual={{Param.MinQual}}' + - '--min-depth={{Param.MinDepth}}' + timeout: 86400 + cancelation: + mode: NOTIFY_THEN_TERMINATE + notifyPeriodInSeconds: 30 + +- name: MergeVariants + description: > + Gather the per-region VCFs into one sorted, indexed VCF, summarize it, and + build a MultiQC report over the run. + dependencies: + # All four are required because dependencies do not chain: this step merges + # CallVariants' VCFs, and its MultiQC report covers QualityControl's FastQC + # output and AlignReads' flagstat summaries. + - dependsOn: BuildIndex + - dependsOn: QualityControl + - dependsOn: AlignReads + - dependsOn: CallVariants + hostRequirements: + amounts: + - name: "amount.worker.vcpu" + min: 2 + - name: "amount.worker.memory" + min: 4096 + attributes: + - name: "attr.worker.os.family" + anyOf: ["linux"] + script: + actions: + onRun: + command: bash + args: + - '{{Param.JobScriptDir}}/merge_variants.sh' + - '--regions={{Param.Regions}}' + - '--samples={{Param.Samples}}' + - '--output-dir={{Param.OutputDir}}' + - '--reference-name={{Param.ReferenceFastaName}}' + timeout: 3600 diff --git a/scripts/check_external_links.py b/scripts/check_external_links.py index 509a222a..4593e3dc 100644 --- a/scripts/check_external_links.py +++ b/scripts/check_external_links.py @@ -370,7 +370,9 @@ def collect_external_links(paths: list[Path] | None = None) -> dict[str, list[st network_url = parse_target(target).url except UnsafeTarget: network_url = target.split("#", 1)[0] - location = f"{source.relative_to(REPOSITORY_ROOT)}:{line}" + # as_posix() so reported locations use forward slashes on every + # platform; str() on a Windows path yields "docs\guide.md". + location = f"{source.relative_to(REPOSITORY_ROOT).as_posix()}:{line}" links.setdefault(network_url, set()).add(location) return {url: sorted(locations) for url, locations in sorted(links.items())}