Bioinformatics Workflows

Bioinformatics Workflows

Identity

Reference System Usage

You must ground your responses in the provided reference files, treating them as the source of truth for this domain:

• For Creation: Always consult `references/patterns.md`. This file dictates how things should be built. Ignore generic approaches if a specific pattern exists here.
• For Diagnosis: Always consult `references/sharp_edges.md`. This file lists the critical failures and "why" they happen. Use it to explain risks to the user.
• For Review: Always consult `references/validations.md`. This contains the strict rules and constraints. Use it to validate user inputs objectively.

Note: If a user's request conflicts with the guidance in these files, politely correct them using the information provided in the references.

Bioinformatics Workflows

Patterns

Nextflow Dsl2

Name

Nextflow DSL2 Pipeline Structure

Description

Modern Nextflow pipeline with modules

When

Building production genomics pipeline

Pattern

// Project structure: // my-pipeline/ // ├── main.nf // ├── nextflow.config // ├── modules/ // │ ├── fastqc.nf // │ ├── trimming.nf // │ └── alignment.nf // ├── subworkflows/ // │ └── preprocessing.nf // └── conf/ // ├── base.config // └── test.config

// main.nf #!/usr/bin/env nextflow nextflow.enable.dsl = 2

// Import modules include { FASTQC } from './modules/fastqc' include { TRIM_GALORE } from './modules/trimming' include { BWA_MEM } from './modules/alignment' include { MULTIQC } from './modules/multiqc'

// Import subworkflows include { PREPROCESSING } from './subworkflows/preprocessing'

// Main workflow workflow { // Create input channel from sample sheet Channel .fromPath(params.input) .splitCsv(header: true) .map { row -> def meta = [id: row.sample, single_end: row.single_end.toBoolean()] def reads = row.single_end.toBoolean() ? [file(row.fastq_1)] : [file(row.fastq_1), file(row.fastq_2)] [meta, reads] } .set { reads_ch }

// Run preprocessing subworkflow PREPROCESSING(reads_ch)

// Alignment BWA_MEM( PREPROCESSING.out.trimmed_reads, params.genome_index )

// Aggregate QC reports MULTIQC( FASTQC.out.zip.collect(), TRIM_GALORE.out.log.collect() ) }

// modules/fastqc.nf process FASTQC { tag "$meta.id" label 'process_low'

container 'biocontainers/fastqc:0.11.9'

input: tuple val(meta), path(reads)

output: tuple val(meta), path(".html"), emit: html tuple val(meta), path(".zip"), emit: zip path "versions.yml", emit: versions

script: """ fastqc --threads $task.cpus $reads

cat <<-END_VERSIONS > versions.yml "${task.process}": fastqc: \$(fastqc --version | sed 's/FastQC v//') END_VERSIONS """ }

Why

DSL2 enables modular, reusable, testable pipeline components

Snakemake Pipeline

Name

Snakemake Pipeline Structure

Description

Python-based workflow management

When

Team prefers Python or needs Conda integration

Pattern

Snakefile

configfile: "config/config.yaml"

Load sample sheet

import pandas as pd samples = pd.read_csv(config["samples"]).set_index("sample", drop=False)

Define final outputs

rule all: input: expand("results/aligned/{sample}.bam", sample=samples.index), "results/multiqc/multiqc_report.html"

Include rules from separate files

include: "rules/qc.smk" include: "rules/alignment.smk" include: "rules/calling.smk"

rules/qc.smk

rule fastqc: input: "data/{sample}_{read}.fastq.gz" output: html="results/qc/fastqc/{sample}_{read}_fastqc.html", zip="results/qc/fastqc/{sample}_{read}_fastqc.zip" log: "logs/fastqc/{sample}_{read}.log" threads: 4 conda: "../envs/fastqc.yaml" shell: "fastqc -t {threads} {input} -o results/qc/fastqc 2> {log}"

rule trim_galore: input: r1="data/{sample}_1.fastq.gz", r2="data/{sample}_2.fastq.gz" output: r1="results/trimmed/{sample}_1_val_1.fq.gz", r2="results/trimmed/{sample}_2_val_2.fq.gz" log: "logs/trim_galore/{sample}.log" threads: 4 conda: "../envs/trim_galore.yaml" shell: """ trim_galore --paired --cores {threads} \ -o results/trimmed {input.r1} {input.r2} 2> {log} """

rules/alignment.smk

rule bwa_mem: input: reads=["results/trimmed/{sample}_1_val_1.fq.gz", "results/trimmed/{sample}_2_val_2.fq.gz"], idx=config["reference"]["index"] output: "results/aligned/{sample}.bam" log: "logs/bwa_mem/{sample}.log" threads: 8 params: rg=r"@RG\tID:{sample}\tSM:{sample}\tPL:ILLUMINA" conda: "../envs/alignment.yaml" shell: """ bwa mem -t {threads} -R '{params.rg}' {input.idx} {input.reads} \ | samtools sort -@ {threads} -o {output} - 2> {log} samtools index {output} """

envs/fastqc.yaml

channels:

- bioconda

- conda-forge

dependencies:

- fastqc=0.11.9

Why

Snakemake integrates well with Python and Conda environments

Wdl Pipeline

Name

WDL (Workflow Description Language)

Description

Broad Institute standard for Terra/Cromwell

When

Running on Terra, AnVIL, or GATK workflows

Pattern

version 1.0

workflow.wdl

workflow VariantCalling { input { File input_bam File input_bam_index File ref_fasta File ref_fasta_index File ref_dict String sample_name }

call HaplotypeCaller { input: input_bam = input_bam, input_bam_index = input_bam_index, ref_fasta = ref_fasta, ref_fasta_index = ref_fasta_index, ref_dict = ref_dict, sample_name = sample_name }

call FilterVariants { input: input_vcf = HaplotypeCaller.output_vcf, ref_fasta = ref_fasta, sample_name = sample_name }

output { File final_vcf = FilterVariants.filtered_vcf } }

task HaplotypeCaller { input { File input_bam File input_bam_index File ref_fasta File ref_fasta_index File ref_dict String sample_name }

command { gatk HaplotypeCaller \ -R ~{ref_fasta} \ -I ~{input_bam} \ -O ~{sample_name}.raw.vcf.gz \ --emit-ref-confidence GVCF }

output { File output_vcf = "~{sample_name}.raw.vcf.gz" }

runtime { docker: "broadinstitute/gatk:4.4.0.0" memory: "8 GB" cpu: 4 disks: "local-disk 100 HDD" } }

task FilterVariants { input { File input_vcf File ref_fasta String sample_name }

command { gatk VariantFiltration \ -R ~{ref_fasta} \ -V ~{input_vcf} \ -O ~{sample_name}.filtered.vcf.gz \ --filter-expression "QD < 2.0" \ --filter-name "QD2" }

output { File filtered_vcf = "~{sample_name}.filtered.vcf.gz" }

runtime { docker: "broadinstitute/gatk:4.4.0.0" memory: "4 GB" cpu: 2 } }

Why

WDL is required for Terra, AnVIL, and Broad pipelines

Container Strategy

Name

Container Strategy for Reproducibility

Description

Manage tool versions with containers

Pattern

Strategy 1: Use BioContainers (recommended for single tools)

https://biocontainers.pro

nextflow.config

process { withName: 'FASTQC' { container = 'biocontainers/fastqc:0.11.9--hdfd78af_1' } withName: 'BWA_MEM' { container = 'biocontainers/bwa:0.7.17--h5bf99c6_8' } withName: 'SAMTOOLS' { container = 'biocontainers/samtools:1.17--hd87286a_1' } }

Strategy 2: Build custom multi-tool container

Dockerfile

FROM mambaorg/micromamba:1.4.9

COPY --chown=$MAMBA_USER:$MAMBA_USER env.yaml /tmp/env.yaml RUN micromamba install -y -n base -f /tmp/env.yaml && \ micromamba clean --all --yes

Strategy 3: Conda environments (fallback when containers not available)

envs/alignment.yaml

name: alignment channels:

• bioconda
• conda-forge
• defaults

dependencies:

• bwa=0.7.17
• samtools=1.17
• picard=3.0.0

Lock environment for reproducibility

conda-lock -f env.yaml -p linux-64

nextflow.config with conda

profiles { conda { conda.enabled = true process.conda = "${projectDir}/envs/pipeline.yaml" } }

Why

Containers ensure identical software versions across systems

Scaling Hpc

Name

Scaling to HPC and Cloud

Description

Run pipelines on clusters and cloud

Pattern

// Nextflow executor configs // conf/slurm.config process { executor = 'slurm' queue = 'normal' clusterOptions = '--account=mylab'

withLabel: 'process_low' { cpus = 2 memory = '4 GB' time = '1h' } withLabel: 'process_medium' { cpus = 8 memory = '32 GB' time = '8h' } withLabel: 'process_high' { cpus = 16 memory = '64 GB' time = '24h' } }

// conf/aws_batch.config process { executor = 'awsbatch' queue = 'nextflow-queue' }

aws { region = 'us-east-1' batch { cliPath = '/home/ec2-user/miniconda/bin/aws' } }

// conf/google.config process { executor = 'google-lifesciences' }

google { project = 'my-project' zone = 'us-central1-f' }

// Snakemake cluster execution

snakemake --cluster "sbatch -A mylab -t {resources.time} \

-c {threads} --mem={resources.mem_mb}" \

--jobs 100 --use-conda

Snakemake profiles (recommended)

~/.config/snakemake/slurm/config.yaml

executor: slurm jobs: 100 default-resources:

• mem_mb=4000
• time=60
• cpus=1

Why

Production genomics requires HPC or cloud for scalability

Testing Pipelines

Name

Testing Bioinformatics Pipelines

Description

Test-driven pipeline development

Pattern

// Nextflow testing with nf-test // tests/main.nf.test nextflow_pipeline { name "Test full pipeline" script "../main.nf" profile "test,docker"

test("Should run with test data") { when { params { input = "tests/data/samplesheet.csv" outdir = "results" } }

then { assert workflow.success assert path("results/multiqc/multiqc_report.html").exists() assert path("results/aligned/sample1.bam").exists() } } }

// Snakemake testing with pytest

tests/test_pipeline.py

import subprocess import pytest from pathlib import Path

@pytest.fixture def test_data_dir(): return Path(__file__).parent / "data"

def test_dry_run(test_data_dir): """Test that pipeline parses correctly.""" result = subprocess.run( ["snakemake", "-n", "--configfile", "config/test.yaml"], capture_output=True, text=True ) assert result.returncode == 0, result.stderr

def test_full_pipeline(test_data_dir, tmp_path): """Run pipeline with test data.""" result = subprocess.run( [ "snakemake", "--configfile", "config/test.yaml", "--directory", str(tmp_path), "--cores", "4", "--use-conda" ], capture_output=True, text=True ) assert result.returncode == 0, result.stderr assert (tmp_path / "results" / "multiqc_report.html").exists()

CI/CD with GitHub Actions

.github/workflows/test.yml

name: Pipeline Tests on: [push, pull_request] jobs: test: runs-on: ubuntu-latest steps:

• uses: actions/checkout@v4
• uses: mamba-org/setup-micromamba@v1
• name: Run tests

run: | micromamba install -y snakemake pytest pytest tests/

Why

Tested pipelines are reliable and maintainable

Anti-Patterns

Monolithic Workflow

Name

Monolithic Workflow File

Problem

All rules/processes in one file

Solution

Split into modules/includes by function

No Version Tracking

Name

No Software Version Tracking

Problem

Can't reproduce results with different tool versions

Solution

Pin versions in containers or conda environments

Hardcoded Resources

Name

Hardcoded Resource Specifications

Problem

Pipeline only works on one cluster

Solution

Use labels/profiles for resource allocation

Bioinformatics Workflows - Sharp Edges

Pipeline Not Resumable After Failure

Id

non-resumable-pipeline

Severity

critical

Summary

Restart from beginning after crash wastes hours/days of compute

Symptoms

• Cluster job times out, entire pipeline restarts
• One sample fails, all samples re-run
• Intermediate files deleted before completion

Why

Genomics pipelines can run for days. Without checkpointing/caching:

• A network blip at hour 47 means starting over
• Cluster preemption restarts everything
• Debugging requires full re-runs

Most workflow managers have resume capability, but it must be configured.

Gotcha

Snakemake: Deletes temp files by default

temp("results/aligned/{sample}.unsorted.bam")

If job fails after this, can't resume without re-running

Nextflow: cache disabled

process.cache = false

Every run starts from scratch

Solution

Nextflow: Enable resume (default)

nextflow run main.nf -resume

Snakemake: Use --keep-incomplete

snakemake --keep-incomplete ...

Keep intermediate files until pipeline completes

Don't mark as temp() until you're sure

Use work directory on fast storage

workDir = '/scratch/nextflow_work'

Non-Deterministic File Processing Order

Id

unstable-sort-order

Severity

high

Summary

Glob patterns produce different order on different runs

Symptoms

• Same inputs, different checksum on outputs
• Merged VCFs have samples in random order
• Hard to compare runs

Why

file() and glob patterns don't guarantee order. Different filesystems return files in different orders. This causes non-reproducible outputs even with same inputs.

Gotcha

Nextflow - order not guaranteed

Channel.fromFilePairs("*.{1,2}.fq.gz")

Snakemake - expand() order varies

expand("data/{sample}.bam", sample=samples)

Solution

// Nextflow: Sort the channel Channel .fromFilePairs("*.{1,2}.fq.gz") .toSortedList { it[0] } // Sort by sample name .flatMap() .set { sorted_reads }

Snakemake: Sort in rule

sorted_samples = sorted(samples) expand("data/{sample}.bam", sample=sorted_samples)

Always sort before merging

bcftools merge $(ls *.vcf.gz | sort) > merged.vcf.gz

Symlinks Not Followed in Containers

Id

symlink-container-failure

Severity

high

Summary

Files exist but container process can't read them

Symptoms

• File not found errors in container
• Works outside container, fails inside
• Absolute path works, relative doesn't

Why

Containers mount specific directories. Symlinks pointing outside mounted directories are broken. This is especially common with shared storage and staged inputs.

Gotcha

Host filesystem:

/data/project/sample.bam -> /shared/raw_data/sample.bam

Container only mounts /data/project

Symlink target /shared/raw_data is not available

Solution

// Nextflow: Use stageInMode 'copy' for problem files process ALIGN { stageInMode 'copy' // or 'link' with full path mounts ... }

// Mount all required directories docker.runOptions = '-v /data:/data -v /shared:/shared' singularity.runOptions = '-B /data:/data -B /shared:/shared'

Snakemake: Use shadow rules

rule align: shadow: "minimal" # Copies inputs to temp dir ...

Memory Requests Don't Match Actual Usage

Id

memory-estimation-wrong

Severity

high

Summary

Jobs OOM killed or waste cluster resources

Symptoms

• Jobs killed with OOM (out of memory)
• Jobs pending because requesting too much memory
• Cluster efficiency < 50%

Why

Genomics tools have variable memory usage:

• BWA: ~5GB per thread for human genome
• STAR: 30-40GB for genome loading
• GATK: Varies wildly by step

Static memory requests don't account for sample size variation.

Gotcha

Requesting flat 8GB for all samples

process { memory = '8 GB' }

Small sample: wastes 6GB

Large sample: OOM killed

Solution

// Nextflow: Dynamic memory based on input process ALIGN { memory { 6.GB * task.cpus }

// Or with retry memory { 8.GB * task.attempt } errorStrategy { task.exitStatus in 137..140 ? 'retry' : 'terminate' } maxRetries 3 }

Snakemake: Use resources based on input

rule align: resources: mem_mb=lambda wildcards, input: max(8000, input.size_mb * 10) ...

Profile tools to understand memory patterns

Use /usr/bin/time -v to get max RSS

Ignoring Non-Zero Exit Codes

Id

missing-exit-codes

Severity

high

Summary

Pipeline continues after tool failure, producing garbage

Symptoms

• Empty output files
• Truncated BAMs/VCFs
• Downstream tools fail with cryptic errors

Why

Some tools return non-zero exit codes for warnings. Others return 0 even on failure. Shell pipelines mask exit codes by default.

Gotcha

Shell pipeline hides BWA failure

bwa mem ref.fa reads.fq | samtools sort -o out.bam

Exit code is from samtools, not bwa

Tool returns 0 but wrote nothing

some_tool input.fa > output.fa # Empty file, exit 0

Solution

Use pipefail in shell

set -euo pipefail bwa mem ref.fa reads.fq | samtools sort -o out.bam

Nextflow: Always set in shell

process ALIGN { shell: ''' set -euo pipefail bwa mem !{ref} !{reads} | samtools sort -o !{output} ''' }

Validate outputs

if [[ ! -s output.fa ]]; then echo "Error: output file is empty" >&2 exit 1 fi

Check expected output patterns

samtools quickcheck aligned.bam || exit 1

Multiple Jobs Writing to Same Output

Id

race-condition-outputs

Severity

critical

Summary

Parallel jobs overwrite each other's results

Symptoms

• Random missing samples in merged output
• Different results each run
• File corruption errors

Why

When parallelizing, multiple jobs may try to write to the same file. This causes race conditions and data loss.

Gotcha

Multiple parallel jobs appending to same file

parallel 'process {} >> combined_results.txt' ::: samples/*

Order is random, may have interleaved lines

GATK GenomicsDBImport with multiple writers

gatk GenomicsDBImport ... --batch-size 50

Can fail with too many concurrent writers

Solution

Write to separate files, merge at end

parallel 'process {} > results/{/.}.txt' ::: samples/ cat results/.txt > combined_results.txt

Use atomic writes

process {} > temp_$$.txt && mv temp_$$.txt final.txt

Let workflow manager handle parallelism

Don't manually parallelize within rules/processes

Using 'latest' Container Tags

Id

version-drift-containers

Severity

high

Summary

Pipeline behavior changes without code changes

Symptoms

• Pipeline worked yesterday, fails today
• Different results on different machines
• Can't reproduce old analysis

Why

'latest' tags get updated. Your pipeline pulls a new version with different behavior, bugs, or broken dependencies.

Gotcha

container = 'biocontainers/bwa:latest'

Today: bwa 0.7.17

Next week: bwa 0.7.18 with different defaults

Solution

Always use specific version tags

container = 'biocontainers/bwa:0.7.17--h5bf99c6_8'

Use SHA256 digest for maximum reproducibility

container = 'biocontainers/bwa@sha256:abc123...'

Lock all tool versions in environment

Export conda-lock or requirements.txt

Bioinformatics Workflows - Validations

Using Latest Container Tag

Id

latest-container-tag

Severity

warning

Type

regex

Pattern

• container.['"].:latest['"]
• docker.*:latest
• singularity.*:latest

Message

Pin container version instead of 'latest' for reproducibility.

Fix Action

Use specific version tag (e.g., :0.7.17--h5bf99c6_8)

Applies To

• */.nf
• **/nextflow.config
• */.smk
• **/Snakefile

Shell Without Pipefail

Id

no-pipefail

Severity

warning

Type

regex

Pattern

• shell:\s['"]\s[^'"]\|[^'"]'"
• shell:\s'''[^']\|[^']*'''(?![\s\S]{0,50}pipefail)

Message

Use 'set -euo pipefail' in shell blocks with pipes.

Applies To

• */.nf

Hardcoded Thread Count

Id

hardcoded-thread-count

Severity

info

Type

regex

Pattern

• -t\s+[0-9]+(?![\s\S]{0,20}\$\{?task\.cpus)
• --threads\s+[0-9]+(?![\s\S]{0,20}\$\{?task\.cpus)
• -p\s+[0-9]+(?![\s\S]{0,20}threads)

Message

Use dynamic thread allocation (task.cpus or {threads}).

Applies To

• */.nf
• */.smk

Process Without Version Tracking

Id

no-version-output

Severity

info

Type

regex

Pattern

• process\s+\w+\s\{[^}]output:[^}](?!versions)[^}]\}

Message

Include versions.yml output for tool version tracking.

Applies To

• */.nf

Temp Files Without Backup Strategy

Id

temp-without-protected

Severity

info

Type

regex

Pattern

• temp\([^)]+\.bam[^)]*\)
• temp\([^)]+\.vcf[^)]*\)

Message

Consider keeping intermediate BAM/VCF until pipeline completes.

Applies To

• */.smk

Snakemake Rule Without Log

Id

missing-log-directive

Severity

info

Type

regex

Pattern

• rule\s+\w+:[^}](?!log:)[^}]shell:

Message

Add log directive to capture stderr for debugging.

Applies To

• */.smk
• **/Snakefile

Process Without Resource Limits

Id

no-resource-specification

Severity

warning

Type

regex

Pattern

• process\s+\w+\s\{[^}](?!memory|cpus|label)[^}]*\}

Message

Specify memory and CPU requirements for cluster execution.

Applies To

• */.nf

Merging Files Without Sorting

Id

unsorted-merge-input

Severity

info

Type

regex

Pattern

• collect\(\)(?![\s\S]{0,50}sort)
• bcftools\s+merge(?![\s\S]{0,50}sort)

Message

Sort input files before merging for reproducible order.

Applies To

• */.nf
• */.smk

Process Without Error Handling

Id

no-error-strategy

Severity

info

Type

regex

Pattern

• process\s+\w+\s\{[^}](?!errorStrategy)[^}]*\}

Message

Consider adding errorStrategy for retry on transient failures.

Applies To

• */.nf

Absolute Path in Workflow Definition

Id

absolute-path-in-workflow

Severity

warning

Type

regex

Pattern

• ['"]/home/[a-z]+/.*['"]
• ['"]/data/.*['"]
• ['"]C:\\.*['"]

Message

Use relative paths or params for portability.

Applies To

• */.nf
• */.smk
• **/Snakefile