Foundations bash scripting

Bash Scripting for Bioinformatics

When to Use

• Automating the same command (FastQC, samtools, BLAST, alignment) over dozens–hundreds of samples.
• Writing the shell glue that Snakemake/Nextflow rules or a cron job calls under the hood.
• Building a CLI wrapper script that validates input files/arguments before running an expensive tool.
• Debugging a pipeline step that "worked once" but silently produced empty/wrong output (missing set -euo pipefail, unquoted variables, unmatched glob).
• Generating a sample sheet or summary report from a directory of FASTQ/BAM/VCF files with awk/sed/cut.

Version Compatibility

• GNU Bash ≥ 4.4 (associative arrays, ${var,,}); Bash 5.x on most modern Linux distros. macOS ships Bash 3.2 — install Bash 5 via Homebrew if targeting that.
• Coreutils/awk/sed/grep as found on any standard Linux distro (GNU variants; BSD sed/awk on macOS differ slightly, e.g. sed -i '').

Prerequisites

• Comfortable with the Linux command line (foundations-linux-fundamentals): paths, pipes, redirection, grep/cut/sort.
• No extra packages to install — Bash and coreutils are present on any POSIX system. Tool-specific commands shown below (fastqc, samtools, bcftools, blastn) assume those tools are separately installed and on $PATH.

Goal: produce a script that never silently continues after a failure and never breaks on filenames containing spaces. Approach: start every script with strict-mode + a logger, validate inputs before doing real work, then loop.

#!/bin/bash
set -euo pipefail
# -e: exit on any non-zero command; -u: error on unset variables (catches typos);
# -o pipefail: a failing command inside a pipe fails the whole pipe.

THREADS=8
INPUT_DIR="${1:?Usage: $0 <input_dir> <output_dir>}"
OUTPUT_DIR="${2:?Usage: $0 <input_dir> <output_dir>}"

log() { echo "[$(date '+%H:%M:%S')] $1" >&2; }

require_tool() {
    # Fail fast if a needed binary isn't on PATH.
    local tool="$1"
    command -v "$tool" &>/dev/null || { log "ERROR: $tool not found"; exit 1; }
}

[[ -d "$INPUT_DIR" ]] || { log "ERROR: not a directory: $INPUT_DIR"; exit 1; }
mkdir -p "$OUTPUT_DIR"
require_tool fastqc
log "Starting..."

Goal: batch-process every FASTQ pair in a directory without breaking on an empty glob or a missing mate. Approach: glob + guard, derive R2 from R1 with pattern substitution, skip (don't crash) on missing pairs.

#!/bin/bash
set -euo pipefail

INPUT_DIR="${1:-.}"
OUTPUT_DIR="${2:-qc_results}"
mkdir -p "$OUTPUT_DIR"

count=0
for r1 in "${INPUT_DIR}"/*_R1.fastq.gz; do
    [[ -f "$r1" ]] || { echo "No *_R1.fastq.gz files found in $INPUT_DIR" >&2; exit 1; }

    r2="${r1/_R1/_R2}"                       # pattern substitution: _R1 -> _R2
    sample=$(basename "$r1" _R1.fastq.gz)

    if [[ ! -f "$r2" ]]; then
        echo "WARNING: missing R2 for $sample, skipping" >&2
        continue
    fi

    echo "[$((++count))] fastqc -t 4 -o \"$OUTPUT_DIR\" \"$r1\" \"$r2\""
    # fastqc -t 4 -o "$OUTPUT_DIR" "$r1" "$r2"
done
echo "Processed $count pairs"

Goal: read a tab-separated sample sheet and summarize gene counts, without touching Python. Approach: while IFS=$'\t' read -r ... for line-oriented parsing; awk -F'\t' for column math.

#!/bin/bash
set -euo pipefail

# Sample sheet: sample_id<TAB>condition
while IFS=$'\t' read -r sample_id condition; do
    echo "Sample: $sample_id | Condition: $condition"
done < sample_sheet.tsv

# Genes with mean count across 3 columns > 1000, from a TSV with a header row
awk -F'\t' 'NR>1 { avg=($2+$3+$4)/3; if (avg>1000) print $1, "avg="int(avg) }' gene_counts.tsv

Case statement for routing by file type:

case "$file" in
    *.fastq.gz | *.fq.gz) fastqc "$file" ;;
    *.bam)                samtools flagstat "$file" ;;
    *.vcf | *.vcf.gz)     bcftools stats "$file" ;;
    *.fasta | *.fa)       grep -c "^>" "$file" ;;
    *) echo "Unknown: $file" >&2; exit 1 ;;
esac

Cleanup on exit (temp dirs, partial output) with trap:

TMPDIR=$(mktemp -d)
cleanup() { rm -rf "$TMPDIR"; }
trap cleanup EXIT     # runs on normal exit AND on error (with set -e)

Comparison Table

Goal	Command
Default value	${var:-default}
Require arg	${1:?Usage: ...}
Strip suffix	${var%.fastq.gz}
Strip extension via basename	$(basename "$f" .gz)
Dir of file	$(dirname "$f")
Count lines	wc -l < "$file"
Redirect stderr	cmd 2>/dev/null
N-way parallel	find . -name '*.gz' | xargs -P4 -I{} fastqc {}

Pitfalls

• set -euo pipefail is non-negotiable: without it, a failed command inside a pipe (cat missing | wc -l) reports the exit code of wc, not cat, and the script continues on garbage.
• No spaces around =: var=value is correct; var = value runs a command named var.
• Always double-quote variables: "$var" not $var — a filename with a space breaks unquoted expansion into multiple words/arguments.
• $() not backticks: backticks cannot be nested and are harder to read.
• local in functions: undeclared variables leak into global scope and can clobber an outer variable with the same name.
• Glob matching zero files: for f in dir/*.gz — when nothing matches, the loop body still runs once with the literal string dir/*.gz; always guard with [[ -f "$f" ]] || { ...; exit 1; } at the top of the loop.
• -u plus optional args: with set -u, referencing $1 when no argument was passed is a hard error — use ${1:-default} or ${1:?message} instead of bare $1.

See Also

• foundations-linux-fundamentals — command-line basics this skill builds on.
• bio-workflow-management-snakemake-workflows — when a script grows past a few pipeline steps, wrap it in Snakemake instead of chaining more Bash.
• bio-workflow-management-nextflow-pipelines — alternative workflow manager for multi-step, multi-sample pipelines.
• foundations-statistics-python — once data leaves the shell (counts, tables), continue analysis in Python/pandas.