Bioinformatics workflows cicd

Bioinformatics Workflows and CI/CD

When to Use

• Turning a tangle of shell scripts into a resumable, parallel DAG pipeline (Snakemake or Nextflow)
• Scaling a pipeline across samples on HPC/cloud (SLURM, AWS Batch, Google Batch)
• Running or customizing an nf-core community pipeline (rnaseq, sarek, atacseq, ...)
• Containerizing pipeline steps with conda environments or Docker/Singularity
• Adding pytest unit tests and a GitHub Actions CI workflow to a bioinformatics repo

Version Compatibility

• Snakemake ≥8.0 (new --executor plugin API replaces --cluster/--cluster-config)
• Nextflow ≥23.10, DSL2 only (DSL1 is deprecated); nf-core pipelines pinned by -r <version>
• pytest ≥7.0, pytest-cov ≥4.0
• GitHub Actions: actions/checkout@v4, actions/setup-python@v5

Prerequisites

conda create -n snake -c conda-forge -c bioconda snakemake
conda install -c bioconda nextflow   # or: curl -s https://get.nextflow.io | bash
pip install pytest pytest-cov nf-core

Prior concepts: shell scripting basics, YAML/config files, Docker/Singularity fundamentals.

Goal: replace a linear shell script with a DAG that resumes after failure, parallelizes independent samples, and pins resources per step.

Approach: define one Snakemake rule per pipeline step with explicit input/output/threads/resources, drive the whole run from a config.yaml, and let rule all declare final targets (Snakemake resolves the DAG backward from there).

## Snakefile
configfile: "config/config.yaml"

SAMPLES = config["samples"]
REF = config["reference"]

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

rule fastqc:
    input: "data/{sample}_R1.fastq.gz"
    output: "results/qc/{sample}_fastqc.html"
    conda: "envs/qc.yaml"
    log: "logs/fastqc/{sample}.log"
    shell: "fastqc {input} --outdir results/qc/ 2> {log}"

rule bwa_mem:
    input:
        r1 = "trimmed/{sample}_R1.fq.gz",
        r2 = "trimmed/{sample}_R2.fq.gz",
        ref = REF
    output: bam = "aligned/{sample}.bam"
    threads: 8
    resources: mem_mb = 16000
    benchmark: "benchmarks/bwa_mem/{sample}.tsv"
    log: "logs/bwa_mem/{sample}.log"
    shell:
        "(bwa mem -t {threads} {input.ref} {input.r1} {input.r2} "
        "| samtools sort -o {output.bam}) 2> {log}"

snakemake -n --configfile config/config.yaml           # dry run
snakemake --cores 8 --use-conda                        # local run, per-rule conda envs
snakemake --dag | dot -Tpng > dag.png                   # visualize DAG
snakemake --executor slurm --jobs 100 --use-conda \
    --default-resources slurm_account=mylab            # HPC (Snakemake >=8 executor API)
snakemake --forcerun bwa_mem --cores 8                  # force re-run one rule

Goal: express the same pipeline as an nf-core-style DSL2 workflow, or run an existing nf-core pipeline instead of writing rules from scratch.

Approach: wrap each tool in a process with a pinned container, wire processes together with channels in a workflow {} block, and prefer nextflow run nf-core/<pipeline> over reinventing common analyses (rnaseq, sarek, atacseq, chipseq...).

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

params.reads  = "data/*_{R1,R2}.fastq.gz"
params.genome = "ref/hg38.fa"
params.outdir = "results"

process BWA_MEM {
    tag "$sample_id"
    publishDir "${params.outdir}/aligned", mode: 'copy'
    container "quay.io/biocontainers/bwa:0.7.17--h5bf99c6_8"
    cpus 8
    memory '16 GB'

    input:
    tuple val(sample_id), path(reads)
    path genome

    output:
    tuple val(sample_id), path("${sample_id}.bam"), emit: bam

    script:
    def (r1, r2) = reads
    """
    bwa mem -t $task.cpus $genome $r1 $r2 | samtools sort -o ${sample_id}.bam
    """
}

workflow {
    reads_ch = Channel.fromFilePairs(params.reads, checkIfExists: true)
    genome_ch = Channel.fromPath(params.genome)
    BWA_MEM(reads_ch, genome_ch)
}

nextflow run main.nf -with-docker
nextflow run main.nf -resume                            # hash-based checkpoint restart
nf-core download rnaseq --revision 3.14.0 --container singularity
nextflow run nf-core/sarek -profile docker -r 3.4.4 \
    --input samplesheet.csv --genome GATK.GRCh38 \
    --tools mutect2,strelka --outdir results/somatic/

Goal: catch off-by-one and edge-case bugs (0- vs 1-based coordinates, empty sequences, ambiguous bases) before they reach a pipeline run, and gate merges on those tests via CI.

Approach: write small pure functions with docstrings, cover them with pytest fixtures and parametrize, then run the suite on every push with GitHub Actions.

## bio_utils.py
def gc_content(sequence: str) -> float:
    """Return GC percentage (0-100) of a DNA sequence; empty input -> 0.0."""
    if not sequence:
        return 0.0
    seq = sequence.upper()
    gc = seq.count("G") + seq.count("C")
    return (gc / len(seq)) * 100

def find_motif(sequence: str, motif: str) -> list[int]:
    """Return all 1-based, overlapping start positions of motif in sequence."""
    positions = []
    start = 0
    while True:
        pos = sequence.find(motif, start)
        if pos == -1:
            break
        positions.append(pos + 1)
        start = pos + 1
    return positions

## test_bio_utils.py
import pytest
from bio_utils import gc_content, find_motif

@pytest.fixture
def sample_fasta_file(tmp_path):
    p = tmp_path / "test.fasta"
    p.write_text(">seq1\nATGCATGC\n>seq2\nGCGCGCGC\n")
    return p

@pytest.mark.parametrize("seq,expected", [
    ("GGGG", 100.0), ("AAAA", 0.0), ("ATGC", 50.0), ("", 0.0),
])
def test_gc_parametrized(seq, expected):
    assert gc_content(seq) == expected

def test_motif_overlapping():
    assert find_motif("AAAA", "AA") == [1, 2, 3]   # overlapping matches, 1-based

## .github/workflows/tests.yml
on: [push, pull_request]
jobs:
  test:
    runs-on: ubuntu-latest
    strategy:
      matrix: { python-version: ["3.10", "3.11", "3.12"] }
    steps:
    - uses: actions/checkout@v4
    - uses: actions/setup-python@v5
      with: { python-version: "${{ matrix.python-version }}" }
    - run: pip install pytest pytest-cov && pip install -r requirements.txt
    - run: pytest --cov=src --cov-report=xml --cov-report=term-missing

Pitfalls

• Snakemake wildcards: rule all must list concrete final targets — Snakemake resolves the DAG backward from outputs, not forward from inputs.
• Snakemake config: always pass --configfile config.yaml; config["samples"] raises KeyError if it's missing or a key is misspelled.
• Nextflow -resume: requires the same work directory and unchanged process definitions (hash-based caching); use it after any partial failure — rerunning without it reprocesses everything.
• Nextflow DSL2 channels: fromFilePairs expects exactly two files per glob match; an odd file count silently drops that sample instead of erroring.
• Container/tool versions: pin tags everywhere (bwa:0.7.17, not latest) — floating tags break reproducibility silently, often months later.
• pytest coordinates: test 0-based vs 1-based boundaries explicitly — off-by-one is the most common bug in bioinformatics test suites.
• CI secrets: never hard-code API keys or tokens; reference them via ${{ secrets.MY_KEY }}.

See Also

• bio-workflow-management-snakemake-workflows — deeper Snakemake patterns (checkpoints, modules)
• bio-workflow-management-nextflow-pipelines — deeper nf-core/Nextflow patterns
• bio-workflow-management-cwl-workflows — CWL as a portable alternative
• clinical-modeling-workflows — ACMG variant classification, Scanpy, docking pipelines