5  The command line as a working environment

After working through this chapter you will be able to

  • navigate a Unix filesystem with absolute and relative paths and read the permission bits on a file
  • connect to a remote server and move data to and from it without retyping passwords or losing work when the connection drops
  • compose grep, sed, awk, sort, uniq and join into pipelines that answer real questions about FASTA, FASTQ, GFF3 and VCF files
  • stream compressed files without decompressing them, and index a VCF for random access
  • diagnose a silently failing pipeline using exit status and set -euo pipefail
  • write a small shell script that takes arguments, validates them and fails loudly
  • submit a job to a SLURM scheduler and request resources that match what the job actually needs

Almost every program a bioinformatician runs is a command-line program. There is no window to click in bwa mem, no menu in bcftools, and no graphical interface on the compute node where a 40-sample alignment will run overnight. That is not an accident of history. A graphical interface is built around one object at a time; the work in this field is built around thousands of objects — reads, genes, variants, samples — and around the composition of small transformations into a pipeline that can be written down, inspected and run again next year on a different dataset.

This chapter treats the shell as a programming environment rather than a list of commands to memorise. The central idea is the stream: a program reads text, transforms it, and writes text another program can read in turn. Once you hold that idea, grep, cut, sort, awk and the rest stop being thirty separate tools and become one toolkit with thirty attachments.

Everything here is worked against the data bundled with this book. You will count reads in data/reads/mgen_sim_R1.fastq.gz, measure sequence lengths in data/proteins/p53-orthologs.fasta, tabulate feature types in the Mycoplasmoides genitalium annotation data/annotation/mycoplasma-genitalium.gff, summarise the clinical classes in the TP53 ClinVar slice data/variants/tp53_clinvar.vcf, and join a protein interaction table to a Gene Ontology annotation file. By the end you will have built a small, reusable command-line tool and run it over every FASTA file in the bundle in parallel.

Each command block is followed by a block holding the real output of that command. All commands are run from the book’s root directory — the one containing data/ — so that every path is relative to a single, predictable place.

5.1 Why the command line

Three properties make the shell the right environment for this work, and none of them is about speed of typing.

The first is composability. A graphical tool does what its author anticipated. A command-line tool does one thing and hands its result to the next, so you can assemble an analysis its authors never imagined. The classic demonstration is a text-processing task, given to a programmer as a week’s work, that Doug McIlroy solved in six lines of shell by piping existing tools together — an exchange that arose from the same body of work as Knuth’s writing on literate programming (Knuth 1984). Here, it means going from a VCF to a table of variant classes without writing a program.

The second is reproducibility. A command is a record of what you did; a sequence of clicks is not. Recording commands in a script, under version control, in a project with a predictable layout is the single highest-return habit in computational research (Wilson et al. 2017; Sandve et al. 2013).

The third is scale. The machine in front of you probably has 8 or 16 GB of memory, and a single whole-genome sequencing run does not fit in it. Command-line tools are built to stream — holding a few lines in memory while terabytes flow past — and the shell is how you drive them. When the data outgrows the laptop, you run the same commands on a server, which has no graphical interface at all.

A shell is a program that reads what you type, works out which program you meant, starts it, and shows you what it wrote. The default on most Linux systems and on macOS is bash or zsh, both of which accept the same core syntax used in this chapter. The shell is also a programming language, with variables, conditionals, loops and functions, which is why a session can be turned into a script with very little editing.

5.2 The shell, the filesystem and permissions

A Unix filesystem is a single tree. Its root is /, and everything — your home directory, an attached hard disk, a shared project volume on a cluster — hangs somewhere beneath it. A path that starts at the root, such as /home/ahmad/projects/pathogens/data, is absolute: it means the same thing from anywhere. A path that does not, such as data/reads, is relative: it is interpreted from your current directory. Two shorthands appear constantly: . is the current directory and .. is its parent.

Relative paths are what make a project portable. If every command in your scripts says data/reads/mgen_sim_R1.fastq.gz, the project can be copied to a cluster and will still work; if one says /Users/ahmad/Desktop/data/reads/mgen_sim_R1.fastq.gz, it works on exactly one machine, for as long as that desktop is tidy.

A dozen commands cover almost all day-to-day navigation and file handling.

Command What it does
pwd Print the absolute path of the current directory
cd DIR Change to DIR; cd alone returns to your home directory
ls, ls -l, ls -lh List files; long form; long form with human-readable sizes
mkdir -p DIR Create a directory, including any missing parents
cp, mv, rm Copy, move or rename, delete
less FILE Page through a file without loading it all (q quits, / searches)
head, tail, wc First lines, last lines, counts of lines, words and bytes
du -sh DIR, df -h . Size of a directory; free space on the filesystem holding it

Begin by looking at what the book ships with.

ls data
MANIFEST.tsv
annotation
expression
networks
proteins
reads
sequences
structures
variants

The long listing adds the information you need when something goes wrong: permissions, owner, size and modification time.

ls -lh data/sequences
total 732K
-rw-r--r-- 1 umarahmad staff 2.6K Oct 15 23:09 TP53-NM_000546.6.fasta
-rw-r--r-- 1 umarahmad staff  37K Oct 15 23:09 TP53-NM_000546.6.gb
-rw-r--r-- 1 umarahmad staff 575K Oct 15 23:09 mycoplasma-genitalium-NC_000908.2.fasta
-rw-r--r-- 1 umarahmad staff  30K Oct 15 23:09 sars-cov-2-NC_045512.2.fasta
-rw-r--r-- 1 umarahmad staff  77K Oct 15 23:09 sars-cov-2-NC_045512.2.gb

The first field is the permission string. Its leading character is the file type — - for an ordinary file, d for a directory, l for a symbolic link — and the next nine are three groups of three: read, write and execute, for the owner, the group and everyone else. So -rw-r--r-- means the owner may read and write and everyone else may only read. For a directory, execute permission means the right to enter it; for a script, the right to run it. chmod +x script.sh adds it, and is the step people forget when their own script reports “permission denied”.

Two more commands answer questions you will ask constantly: how big is this, and what actually is it?

du -sh data
file data/sequences/sars-cov-2-NC_045512.2.fasta data/reads/mgen_sim_R1.fastq.gz
16M data
data/sequences/sars-cov-2-NC_045512.2.fasta: ASCII text
data/reads/mgen_sim_R1.fastq.gz:             gzip compressed data, last modified: Thu Oct 15 20:18:15 2026, original size modulo 2^32 6644844

file inspects content rather than trusting the extension, which matters more than it sounds: a file named reads.fastq that file reports as gzip compressed data is the cause of a surprising share of “the aligner crashed immediately” problems.

The last routine check is integrity. Data that travelled over a network can arrive truncated, and a truncated FASTQ will produce an analysis that runs to completion and is wrong. Every file in this book’s bundle has a recorded SHA-256 checksum in data/MANIFEST.tsv, so the file you have can be compared with the file the author shipped.

sha256sum data/sequences/sars-cov-2-NC_045512.2.fasta
grep 'sars-cov-2-NC_045512.2.fasta' data/MANIFEST.tsv | cut -f3
0891c00cc2d503f58e0c600ce5ee64e9c64b939525c70d5486bc46cbc7cb073e  data/sequences/sars-cov-2-NC_045512.2.fasta
0891c00cc2d503f58e0c600ce5ee64e9c64b939525c70d5486bc46cbc7cb073e

The two strings match, so the file is byte-for-byte what it should be. On macOS the equivalent command is shasum -a 256.

Deleting and overwriting happen silently

rm does not move files to a wastebasket. There is no undo, and on a shared server there may be no backup either. Two habits prevent most disasters: never type rm -rf with a variable in the path until you have run the same command with echo in front of it and read what it would have deleted, and never run rm -rf $DIR/ when $DIR might be empty — rm -rf / is what that becomes.

Redirection with > is equally silent. sort data.txt > data.txt does not sort the file; the shell truncates the target to zero bytes before sort ever opens it, and the data is gone. Write to a new name, check the result, then rename. Adding set -o noclobber to your shell configuration makes > refuse to overwrite an existing file, and >| overrides it when you mean it.

5.3 Working on a remote machine

Most real analysis happens on a machine you do not sit in front of. Four tools cover the whole workflow; what they print depends on the server you connect to, so no output is reproduced here.

You open a shell on the remote machine with ssh, giving a username and a hostname:

ssh ahmad@hpc.example.ac.uk

Typing a password every time is tedious and weaker than the alternative. Generate a key pair once, keep the private half on your laptop and install the public half on the server:

ssh-keygen -t ed25519 -C "ahmad@laptop"
ssh-copy-id ahmad@hpc.example.ac.uk

Moving data is scp for one-off copies and rsync for anything you will repeat; the colon separates the host from the remote path:

scp results/variants.vcf.gz ahmad@hpc.example.ac.uk:/scratch/ahmad/project/
rsync -avz --progress /scratch/ahmad/project/results/ ahmad@hpc.example.ac.uk:backup/results/

rsync transfers only the differences between source and destination, so re-running it after adding three files copies three files rather than three hundred. Note the trailing slash on the source: results/ means “the contents of results”, while results means “the directory itself”. Adding --dry-run shows what would be transferred without transferring it, which is worth doing the first time you use --delete.

The fourth tool solves a problem every remote user meets within a week. If you start a six-hour alignment over ssh and your laptop sleeps, the connection drops and the job dies with it. A terminal multiplexer keeps the session alive on the server, detached from your connection:

tmux new -s alignment
# start the long job, then press Ctrl-b then d to detach
tmux attach -t alignment

screen does the same with different keystrokes (Ctrl-a d to detach, screen -r to reattach) and is installed on some systems where tmux is not. Either way the rule is the same: anything that will run for more than a few minutes starts inside a multiplexer, or inside a scheduler job (Section 5.9).

5.4 Streams, pipes and redirection

Every Unix program starts with three channels open. Standard input (stdin) is where it reads from, standard output (stdout) is where its results go, and standard error (stderr) is where its complaints go. By default stdin is your keyboard and both outputs are your screen, which is why the distinction is invisible until you need it.

Separating results from complaints is the point. A program can write a good table to stdout while warning about a malformed record on stderr, and you can capture the table without the warnings contaminating it.

ls data/sequences/*.fasta data/proteins/missing.fasta > work/found.txt 2> work/errors.txt
cat work/found.txt
cat work/errors.txt
data/sequences/TP53-NM_000546.6.fasta
data/sequences/mycoplasma-genitalium-NC_000908.2.fasta
data/sequences/sars-cov-2-NC_045512.2.fasta
ls: cannot access 'data/proteins/missing.fasta': No such file or directory

The notation is systematic: > redirects stdout to a file, >> appends instead of overwriting, 2> redirects stderr, 2>&1 sends stderr to wherever stdout is currently going, and < file feeds a file to stdin. /dev/null is a file that discards everything written to it, so 2> /dev/null means “silence the warnings” — useful when you know what they are, dangerous when you do not.

The pipe | connects one program’s stdout directly to the next program’s stdin, with no intermediate file. This is the mechanism that turns small tools into an analysis:

grep '^>' data/sequences/*.fasta | tee work/headers.txt | wc -l
3

tee is a useful middle element: it writes its input to a file and passes it on, so you can keep an intermediate result without breaking the pipeline.

Why text, and why line-oriented

FASTA, FASTQ, GFF3, VCF, SAM and BED are all plain text, and all of them are line-oriented: one record, or one part of a record, per line. That is not a historical accident but a design decision, and it is what makes the entire toolkit in this chapter apply to formats invented decades after the tools. A format that is line-oriented text can be counted with wc, filtered with grep, cut into columns with cut, summarised with awk and streamed through a pipe of arbitrary length without ever being loaded into memory.

5.5 The text-processing toolkit

This is the core of the chapter. Each tool below is introduced on a real file and a real question.

Counting and extracting with grep

grep prints lines matching a pattern. In a FASTA file every header begins with >, so counting sequences is counting those lines; -c counts instead of printing.

grep -c '^>' data/proteins/p53-orthologs.fasta
grep '^>' data/proteins/globins.fasta
8
>sp|P69905|HBA_HUMAN Hemoglobin subunit alpha OS=Homo sapiens OX=9606 GN=HBA1 PE=1 SV=2
>sp|P68871|HBB_HUMAN Hemoglobin subunit beta OS=Homo sapiens OX=9606 GN=HBB PE=1 SV=2
>sp|P02144|MYG_HUMAN Myoglobin OS=Homo sapiens OX=9606 GN=MB PE=1 SV=2
>sp|P02185|MYG_PHYMC Myoglobin OS=Physeter macrocephalus OX=9755 GN=MB PE=1 SV=2

The pattern ^> is a regular expression: a small language for describing text patterns. ^ anchors to the start of a line, $ to the end, . matches any character, [A-Z] any character in a range, * zero or more of the preceding item, + one or more (with grep -E), and \t a tab in the tools that accept it. Anchoring matters: grep '>' would also match a > appearing inside a sequence description, and grep 'gene' matches gene, pseudogene and gene_biotype alike. Use -w to match whole words and -i to ignore case.

Two more options earn their place: -v inverts the match, which is how you drop header lines from a VCF or GFF3, and -o prints only the matching part, which turns grep into an extractor:

grep -o -E 'gene=rp[ls][A-Z]' data/annotation/mycoplasma-genitalium.gff | sort -u | head -6
grep -c -i 'ribosomal protein' data/annotation/mycoplasma-genitalium.gff
gene=rplA
gene=rplB
gene=rplC
gene=rplD
gene=rplE
gene=rplF
52

Columns, sorting and counting

Tab-separated formats — GFF3, VCF, BED, the manifest itself — are handled by cut (select columns), sort and uniq -c (count identical adjacent lines). The three together form the most-used idiom in the shell: cut, sort, count.

GFF3 column 3 is the feature type. One pipeline tells you what the annotation actually contains:

grep -v '^#' data/annotation/mycoplasma-genitalium.gff | cut -f3 | sort | uniq -c | sort -k1,1nr
    546 gene
    524 CDS
     42 exon
     36 tRNA
     20 pseudogene
      3 rRNA
      1 RNase_P_RNA
      1 SRP_RNA
      1 region
      1 tmRNA

Read that output as biology, not bookkeeping. M. genitalium has 546 annotated genes, of which 524 are protein-coding, 20 are pseudogenes and the remainder are structural RNAs — a genome stripped to near the minimum a free-living cell can manage. Note also that uniq -c only collapses adjacent duplicates, which is why sort must come first; forgetting that is the most common error with this idiom.

The same pattern applied to the manifest shows how the bundled data is organised:

tail -n +2 data/MANIFEST.tsv | cut -f1 | cut -d/ -f1 | sort | uniq -c
      2 annotation
      2 expression
      1 networks
      4 proteins
      2 reads
      5 sequences
      4 structures
      1 variants

tail -n +2 starts at line 2, discarding the header; cut -d/ -f1 splits on / instead of tab and keeps the first field. Sorting numerically requires care: sort -k1,1nr sorts on field 1 only, numerically (n), in reverse (r). Writing sort -n alone sorts on the whole line, and sort -k1nr sorts from field 1 to the end of the line, which is rarely what you want.

Editing streams with sed

sed applies an edit to every line as it passes. The command worth learning first is substitution, s/pattern/replacement/, with -E for modern regular expression syntax and parenthesised groups recalled as \1, \2. UniProt FASTA headers pack the accession and entry name into a pipe-delimited field; one substitution turns them into a two-column table:

grep '^>' data/proteins/globins.fasta | sed -E 's/^>[a-z]+\|([A-Z0-9]+)\|(\S+).*/\1\t\2/'
P69905  HBA_HUMAN
P68871  HBB_HUMAN
P02144  MYG_HUMAN
P02185  MYG_PHYMC

tr is a cruder relative that translates or deletes individual characters, and it has one classic bioinformatics use — complementing a sequence, which combined with rev gives the reverse complement:

head -2 data/sequences/sars-cov-2-NC_045512.2.fasta | tail -1 | cut -c1-40
head -2 data/sequences/sars-cov-2-NC_045512.2.fasta | tail -1 | cut -c1-40 | tr 'ACGT' 'TGCA' | rev
ATTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCAAC
GTTGGTTGGTTTGTTACCTGGGAAGGTATAAACCTTTAAT

Arithmetic and records with awk

awk reads a file line by line, splits each line into fields, and runs your code on the ones that match a condition. The structure is awk 'CONDITION { ACTION }', fields are $1, $2 and so on with $0 the whole line, NR is the current record number, and -F sets the field separator. It is the right tool as soon as a question involves arithmetic.

FASTA sequence lengths require accumulating across lines until the next header, which is exactly what awk’s state between records is for:

awk '/^>/ {if (name != "") print name "\t" len; name = substr($1, 2); len = 0; next}
     {len += length($0)}
     END {print name "\t" len}' data/proteins/p53-orthologs.fasta
sp|P04637|P53_HUMAN 393
sp|P02340|P53_MOUSE 390
sp|P10361|P53_RAT   391
sp|P41685|P53_FELCA 386
sp|P67939|P53_BOVIN 386
sp|Q9TTA1|P53_TUPBE 393
sp|P79734|P53_DANRE 373
sp|P79820|P53_ORYLA 352

The p53 orthologues range from 352 residues in medaka to 393 in human and tree shrew — a useful sanity check before aligning them, because a sequence far outside that range would signal a fragment or a mis-picked isoform.

For GFF3, awk answers the question “which is the longest gene?” by computing end minus start, with sed tidying the attribute column down to the gene name:

awk -F'\t' '$3 == "gene" {print $5 - $4 + 1 "\t" $9}' data/annotation/mycoplasma-genitalium.gff \
  | sed -E 's/\t.*Name=([^;]+);.*/\t\1/' \
  | sort -k1,1nr | head -5
5418    hmw2
5352    MG_RS02750
4851    MG_RS02355
4356    polC
4335    mgpA

The longest genes in this genome are the cytadherence proteins HMW2 and MgpA and the replicative polymerase PolC. Summary statistics come from the same loop:

awk -F'\t' '$3 == "gene" {len += $5 - $4 + 1; n++}
            END {printf "%d genes, mean length %.1f bp, total %d bp\n", n, len/n, len}' \
    data/annotation/mycoplasma-genitalium.gff
546 genes, mean length 1003.5 bp, total 547911 bp

Gene features span 547,911 bp of a 580,076 bp genome: about 94% of this organism’s DNA is annotated as genes, against roughly 2% of the human genome that is protein-coding.

FASTQ: four lines per read

A FASTQ record is exactly four lines — identifier, sequence, a + separator, and quality characters — which makes NR % 4 the key to the format. Read counts, length distributions and base composition all follow:

zcat data/reads/mgen_sim_R1.fastq.gz | wc -l
echo $(( $(zcat data/reads/mgen_sim_R1.fastq.gz | wc -l) / 4 ))
zcat data/reads/mgen_sim_R1.fastq.gz | awk 'NR % 4 == 2 {print length($0)}' | sort | uniq -c
80000
20000
  20000 150
zcat data/reads/mgen_sim_R1.fastq.gz \
  | awk 'NR % 4 == 2 {total += length($0); gc += gsub(/[GCgc]/, "")}
         END {printf "GC content: %.2f%%\n", 100 * gc / total}'
GC content: 31.70%

The order of the two statements inside the braces is not a matter of taste. gsub returns the number of substitutions it made and modifies $0 in place, deleting every G and C from the line. If total += length($0) came second it would measure the line with those bases already removed, and the reported GC content would be 46.42% — a number that is wrong but entirely plausible, which is the worst kind. The correct figure, 31.70%, agrees with the 31.69% GC of the M. genitalium genome these reads were simulated from, and that agreement is the check worth making: read composition that does not match the reference means contamination, adapter read-through or the wrong reference.

Know when to stop using the shell

paste - - - - converts a four-line FASTQ record into a single tab-separated row, which makes awk filtering on sequence and quality together straightforward:

zcat data/reads/mgen_sim_R1.fastq.gz | paste - - - - \
  | awk -F'\t' '{print $1 "\t" length($2)}' | head -3
@SIM:1:FC1:1:1:0:386803/1   150
@SIM:1:FC1:1:1:1:522982/1   150
@SIM:1:FC1:1:1:2:284660/1   150

For a quick check that is the fastest route. But quality-aware trimming, adapter removal and duplicate marking are solved problems: use Trimmomatic (Bolger, Lohse, and Usadel 2014) or FastQC (Andrews 2010) and aggregate the reports with MultiQC (Ewels et al. 2016). The shell is for asking questions of your data, not for reimplementing tools tested on millions of samples.

Joining two tables

join merges two sorted files on a shared key, which is how you attach annotation to a list of identifiers without opening a spreadsheet. The bundled STRING network lists high-confidence interactions among 30 cell-cycle and DNA-damage genes (Szklarczyk et al. 2022); the UniProt file lists Gene Ontology terms for the reviewed human proteome (Ashburner et al. 2000; Bateman et al. 2022). Both must be sorted on the key column before joining.

tail -n +2 data/networks/p53_string_network.tsv | cut -f3,4 | tr '\t' '\n' | sort -u > work/network_genes.txt
tail -n +2 data/annotation/human_go_annotations.tsv | cut -f2,3 | sort -k1,1 > work/gene_go.tsv
wc -l < work/network_genes.txt
wc -l < work/gene_go.tsv
30
20431
join -t $'\t' work/network_genes.txt work/gene_go.tsv \
  | awk -F'\t' '{print $1 "\t" gsub(/GO:/, "")}' | sort -k2,2nr | head -5
SIRT1   154
AKT1    136
TP53    119
PTEN    91
EGFR    85

The gsub call counts occurrences of GO: and returns that count, so the second column is the number of GO terms annotating each gene. The ordering reflects research attention as much as biology — the most-studied proteins accumulate the most annotations, which is a bias you must allow for in any enrichment analysis.

Two diagnostics belong with every join. join -v 1 prints the keys from file 1 that found no partner, and comparing input and output line counts reveals duplicate keys:

join -t $'\t' work/network_genes.txt work/gene_go.tsv | wc -l
join -t $'\t' work/network_genes.txt work/gene_go.tsv | cut -f1 | uniq -d
31
CDKN2A

Thirty genes went in and 31 rows came out, because CDKN2A has two reviewed UniProt entries — the p16 and p14ARF products of alternative reading frames. A join that silently duplicates rows will silently inflate every count downstream, so always check that the row count is what you expected.

Summarising a VCF

The Variant Call Format packs its annotations into a semicolon-delimited INFO column (Danecek et al. 2011), which is awkward for cut and natural for sed or grep -o. The bundled file is a ClinVar slice over TP53 (Landrum et al. 2017).

grep -c -v '^#' data/variants/tp53_clinvar.vcf
grep -v '^#' data/variants/tp53_clinvar.vcf | cut -f1-5 | head -3
3888
17  7668434 35555   T   G
17  7668440 3051986 T   C
17  7668449 325619  A   C
grep -v '^#' data/variants/tp53_clinvar.vcf \
  | sed -E 's/.*CLNVC=([^;]+).*/\1/' | sort | uniq -c | sort -k1,1nr
   2745 single_nucleotide_variant
    664 Deletion
    239 Duplication
    133 Indel
     52 Insertion
     47 Microsatellite
      8 Inversion

Seventy per cent of the records are single-nucleotide variants, and the remainder are small insertions and deletions — the mutational spectrum expected for a tumour suppressor assayed mostly by exome and panel sequencing.

For anything beyond counting, use a format-aware tool. bcftools query understands the INFO column and takes a format string, which is both shorter and immune to the quoting traps of regular expressions (Danecek et al. 2021):

bcftools query -f '%POS\t%REF\t%ALT\t%INFO/CLNSIG\n' data/variants/tp53_clinvar.vcf | head -5
7668434 T   G   Benign
7668440 T   C   Likely_benign
7668449 A   C   Uncertain_significance
7668477 G   A   Conflicting_classifications_of_pathogenicity
7668539 G   A   Benign/Likely_benign
bcftools query -f '%INFO/MC\n' data/variants/tp53_clinvar.vcf \
  | tr ',' '\n' | cut -d'|' -f2 | sort | uniq -c | sort -k1,1nr | head -6
   1516 missense_variant
    856 non-coding_transcript_variant
    688 intron_variant
    673 5_prime_UTR_variant
    650 frameshift_variant
    513 synonymous_variant

Missense variants dominate, as they must for a gene whose pathogenic mechanism is a dominant-negative DNA-binding domain substitution rather than loss of the protein.

5.6 Compression and streaming

Sequencing data is stored compressed, always. The bundled read file is 2.9 MB on disk and 6.6 MB uncompressed, a ratio that holds roughly for real data too:

gzip -l data/reads/mgen_sim_R1.fastq.gz
         compressed        uncompressed  ratio uncompressed_name
            3010012             6644844  54.7% data/reads/mgen_sim_R1.fastq

The important point is not the disk saving but that you never need to decompress. zcat writes the decompressed stream to stdout, so a compressed file slots into a pipeline where a plain one would, and nothing large is written to disk:

zcat data/reads/mgen_sim_R1.fastq.gz | awk 'NR % 4 == 2' | wc -l
20000

zgrep and zless are the equivalent wrappers for grep and less. Almost every modern bioinformatics tool also reads gzip-compressed input directly, so gunzip before an analysis is usually wasted time and wasted disk.

Plain gzip has one limitation: to reach the middle of the file you must decompress everything before it. The BGZF format solves this by compressing in independent blocks, which is why VCF, BAM and indexed FASTA all use it. bgzip writes it and tabix builds the index (Danecek et al. 2021):

bgzip -c data/variants/tp53_clinvar.vcf > work/tp53_clinvar.vcf.gz
tabix -p vcf work/tp53_clinvar.vcf.gz
ls -lh work/tp53_clinvar.vcf.gz work/tp53_clinvar.vcf.gz.tbi
-rw-r--r-- 1 umarahmad staff 97K Oct 16 01:22 work/tp53_clinvar.vcf.gz
-rw-r--r-- 1 umarahmad staff 151 Oct 16 01:22 work/tp53_clinvar.vcf.gz.tbi

A 1.1 MB VCF becomes 97 KB with a 151-byte index, and any region can now be fetched without touching the rest of the file:

tabix work/tp53_clinvar.vcf.gz 17:7676040-7676050 \
  | awk -F'\t' 'length($4) == 1 && length($5) == 1 {print $1, $2, $4, $5}' | head -5
17 7676040 C A
17 7676040 C G
17 7676040 C T
17 7676041 G A
17 7676041 G C

The awk filter is not decoration. A region query returns every record overlapping the region, and ClinVar contains multi-kilobase deletions whose REF field holds the entire deleted sequence; without the length filter an eleven-base query prints several screens of raw sequence.

5.7 Exit status and defensive scripting

Every command returns an exit status: zero for success, non-zero for failure. The shell stores it in $?, and it is the only reliable way a script can tell whether the previous step worked.

grep -c '^>' data/proteins/no_such_file.fasta; echo "exit status: $?"
grep: data/proteins/no_such_file.fasta: No such file or directory
exit status: 2

Now put that same failing command into a pipeline, which is how it will appear in real work:

grep '^>' data/proteins/no_such_file.fasta | wc -l; echo "exit status: $?"
grep: data/proteins/no_such_file.fasta: No such file or directory
0
exit status: 0

Read that output carefully, because it is the most dangerous behaviour in the shell. grep failed. wc received an empty stream, counted zero lines, and succeeded. The pipeline’s exit status is the status of its last command, so the shell reports success. A script that checks $? sees success and carries on, and you have a results table containing a zero that means “no file” rather than “no sequences”. In a loop over 200 samples, with the error message scrolling past among thousands of other lines, this is how a published count ends up wrong.

The fix is one shell option:

set -o pipefail
grep '^>' data/proteins/no_such_file.fasta | wc -l; echo "exit status: $?"
grep: data/proteins/no_such_file.fasta: No such file or directory
0
exit status: 2

With pipefail set, the pipeline reports the failure of any component. Three options belong at the top of every script you write:

  • set -e — exit immediately if any command fails, instead of carrying on with bad input
  • set -u — treat an unset variable as an error, so a typo in $sampel stops the script rather than expanding to nothing
  • set -o pipefail — as above

They are conventionally written as one line, set -euo pipefail.

An empty variable is a path to somewhere

Without set -u, a mistyped variable expands to the empty string. rm -rf "$OUTDIR"/* with OUTDIR unset becomes rm -rf /*. With set -u the script stops at the unset variable, which is why the three options above are not optional for any script that writes or deletes files.

A script that takes arguments

A script is a file of commands with a line at the top telling the system which interpreter to use. Beyond that, what separates a script you can trust from a sequence of commands in a file is that it checks its inputs and reports failures on stderr with a non-zero exit status. Here is a complete tool that prints the length of every sequence in a FASTA file, with an optional minimum length:

#!/usr/bin/env bash
# Print the length of every sequence in a FASTA file as a two-column table.
set -euo pipefail

usage() {
    echo "usage: $(basename "$0") FASTA [MIN_LENGTH]" >&2
    echo "  writes 'name<TAB>length' to stdout, optionally keeping only" >&2
    echo "  sequences at least MIN_LENGTH residues long" >&2
}

if [[ $# -lt 1 || $# -gt 2 ]]; then
    usage
    exit 2
fi

fasta=$1
min_length=${2:-0}

if [[ ! -r $fasta ]]; then
    echo "$(basename "$0"): cannot read '$fasta'" >&2
    exit 1
fi

awk -v min="$min_length" '
    /^>/ { if (name != "") emit(); name = substr($1, 2); len = 0; next }
          { len += length($0) }
    END  { if (name != "") emit() }
    function emit() { if (len >= min) print name "\t" len }
' "$fasta"

Several details are deliberate. #!/usr/bin/env bash finds bash through the user’s PATH rather than assuming /bin/bash. $# is the number of arguments and $1 the first; ${2:-0} supplies a default when the second is absent. Error messages go to stderr with >&2 so they do not contaminate a redirected table. The exit statuses are conventional: 2 for misuse, 1 for a runtime failure. And "$fasta" is quoted everywhere, because an unquoted variable containing a space becomes two arguments.

Save it as bin/fasta_lengths.sh, make it executable with chmod +x bin/fasta_lengths.sh, and it behaves like any other command:

bin/fasta_lengths.sh data/proteins/globins.fasta
sp|P69905|HBA_HUMAN 142
sp|P68871|HBB_HUMAN 147
sp|P02144|MYG_HUMAN 154
sp|P02185|MYG_PHYMC 154
bin/fasta_lengths.sh data/proteins/p53-orthologs.fasta 390
sp|P04637|P53_HUMAN 393
sp|P02340|P53_MOUSE 390
sp|P10361|P53_RAT   391
sp|Q9TTA1|P53_TUPBE 393

The failure paths matter as much as the success path, and you should test them deliberately:

bin/fasta_lengths.sh
bin/fasta_lengths.sh data/proteins/missing.fasta; echo "exit status: $?"
usage: fasta_lengths.sh FASTA [MIN_LENGTH]
  writes 'name<TAB>length' to stdout, optionally keeping only
  sequences at least MIN_LENGTH residues long
fasta_lengths.sh: cannot read 'data/proteins/missing.fasta'
exit status: 1

Jobs in the background

A command followed by & runs in the background and returns the prompt immediately. jobs lists what is running in the current shell and kill %1 stops job number one:

sleep 20 &
jobs
kill %1
jobs
[1]+  Running                 sleep 20 &
[1]+  Terminated: 15          sleep 20

Ctrl-c interrupts the foreground job, Ctrl-z suspends it, bg resumes it in the background and fg brings it back. For anything that must survive a lost connection, use tmux or a scheduler rather than &. To see how long something took, put time in front of it:

time (zcat data/reads/mgen_sim_R1.fastq.gz | awk 'NR % 4 == 2' | wc -l)
20000

real    0m0.078s
user    0m0.077s
sys 0m0.033s

real is wall-clock time, user is CPU time in your program and sys is CPU time in the kernel. When user exceeds real the job is using multiple cores; when real greatly exceeds user plus sys it is waiting for disk or network, and adding cores will not help.

5.8 Batch work with find, xargs and parallel

find walks a directory tree and selects files by name, type, size or age, and unlike a shell glob it descends through subdirectories:

find data -name '*.fasta'
find data -type f -size +1M -printf '%s\t%p\n' | sort -k1,1nr
data/proteins/globins.fasta
data/proteins/P04637.fasta
data/proteins/p53-orthologs.fasta
data/sequences/sars-cov-2-NC_045512.2.fasta
data/sequences/TP53-NM_000546.6.fasta
data/sequences/mycoplasma-genitalium-NC_000908.2.fasta
5366641 data/annotation/human_go_annotations.tsv
3010341 data/reads/mgen_sim_R2.fastq.gz
3010012 data/reads/mgen_sim_R1.fastq.gz
1126163 data/variants/tp53_clinvar.vcf

xargs turns that list of names into arguments for another command. The pairing of find -print0 with xargs -0 separates names by a null byte rather than whitespace, which is the only safe way to handle filenames containing spaces:

find data -name '*.fasta' -print0 | xargs -0 -n 1 bin/fasta_lengths.sh | sort -k2,2nr | head -4
NC_000908.2 580076
NC_045512.2 29903
NM_000546.6 2512
sp|P04637|P53_HUMAN 393

Your own tool now runs over every FASTA file in the bundle, sorted by sequence length, in one line. When the command needs the filename in the middle rather than at the end, -I {} gives it a placeholder:

find data -name '*.fasta' | xargs -I {} sh -c 'printf "%s\t%s\n" "{}" "$(grep -c "^>" "{}")"'
data/proteins/globins.fasta 4
data/proteins/P04637.fasta  1
data/proteins/p53-orthologs.fasta   8
data/sequences/sars-cov-2-NC_045512.2.fasta 1
data/sequences/TP53-NM_000546.6.fasta   1
data/sequences/mycoplasma-genitalium-NC_000908.2.fasta  1

GNU parallel does the same job across several cores, which is what you want when each task takes minutes rather than milliseconds. ::: introduces the list of inputs, {} is the placeholder, and -j 4 runs four tasks at a time:

parallel -j 4 'printf "%s\t%s\n" {} $(grep -c "^>" {})' ::: data/proteins/*.fasta data/sequences/*.fasta
data/proteins/P04637.fasta  1
data/proteins/globins.fasta 4
data/proteins/p53-orthologs.fasta   8
data/sequences/TP53-NM_000546.6.fasta   1
data/sequences/mycoplasma-genitalium-NC_000908.2.fasta  1
data/sequences/sars-cov-2-NC_045512.2.fasta 1

Note that the order differs from the xargs version: with four workers, whichever finishes first prints first. If downstream code depends on order, sort the output or write one file per input. The first time you run parallel it asks you to cite it; running parallel --citation once records that you have seen the request and silences it.

For anything more complex than “run this command on these files” — a pipeline with dependencies between steps, partial re-runs after a failure, or execution on a cluster — move to a workflow manager such as Snakemake (Köster and Rahmann 2012) or Nextflow (Di Tommaso et al. 2017). They exist precisely because shell loops stop being maintainable at about the point where you want to re-run only the samples that failed.

5.9 Running jobs on a shared cluster

A high-performance computing (HPC) cluster is a set of compute nodes, a shared filesystem and a scheduler that decides which job runs where. You log in to a login node, which is for editing files and submitting work — never for running an aligner. SLURM is the scheduler you are most likely to meet.

Work is submitted as a script whose leading comments are scheduler directives:

#!/usr/bin/env bash
#SBATCH --job-name=align_mgen
#SBATCH --account=pathogen_genomics
#SBATCH --partition=standard
#SBATCH --cpus-per-task=8
#SBATCH --mem=16G
#SBATCH --time=04:00:00
#SBATCH --array=1-24
#SBATCH --output=logs/align_%A_%a.out
#SBATCH --error=logs/align_%A_%a.err

set -euo pipefail

sample=$(sed -n "${SLURM_ARRAY_TASK_ID}p" samples.txt)

bwa mem -t "${SLURM_CPUS_PER_TASK}" reference/mgen.fa \
        "reads/${sample}_R1.fastq.gz" "reads/${sample}_R2.fastq.gz" \
  | samtools sort -@ 2 -o "alignments/${sample}.bam" -
samtools index "alignments/${sample}.bam"

The managing commands are few: sbatch run_alignment.sh submits the script and prints its job identifier; squeue -u $USER lists your queued and running jobs; scancel JOBID cancels one; sacct -j JOBID --format=JobID,Elapsed,MaxRSS,State reports how long a finished job took and how much memory it used.

Several things in that script are worth copying. The job array (--array=1-24) submits 24 independent tasks from one script, each receiving a different SLURM_ARRAY_TASK_ID, which the sed -n "${N}p" line turns into a sample name from a list — this is how per-sample work is parallelised on a cluster. The thread count passed to bwa is ${SLURM_CPUS_PER_TASK} rather than a hard-coded 8, so changing the directive changes the job consistently. Standard output and standard error go to separate files named with both the array job and task identifiers, so that a failure in task 17 can be found. And set -euo pipefail is there for the reason given above: a scheduler reports a job as COMPLETED if its script exits zero, however wrong the contents of the output files are.

Requesting resources well is a social skill as much as a technical one. Ask for what the job needs, measured from a test run with sacct, not what you hope will be enough: a job requesting 128 GB and 48 hours sits in the queue far longer than one requesting 16 GB and four hours, and blocks other people’s work while it runs. Write intermediate files to the cluster’s scratch filesystem rather than your home directory, which is usually small and backed up at everyone’s expense. Never run an alignment on the login node. And install software with Bioconda (Grüning et al. 2018) or the cluster’s module system rather than compiling tools into your home directory and hoping.

Test on one sample before submitting 200

Run the body of the job script interactively on a single sample, in an interactive allocation (srun --pty --cpus-per-task=2 --mem=4G --time=01:00:00 bash), and check the output file. The cost of discovering a typo in an output path after submitting a 200-task array is 200 failed tasks and an apologetic email to the cluster administrator.

5.10 Summary

  • The shell is a programming environment, not a set of commands to memorise. Its organising idea is the stream: programs read text from stdin, write results to stdout and complaints to stderr, and pipes compose them.
  • Relative paths make a project portable; absolute paths tie it to one machine.
  • grep, cut, sort, uniq, sed, awk and join answer most routine questions about FASTA, FASTQ, GFF3 and VCF files, because all of these formats are line-oriented text.
  • sort before uniq -c; check row counts before and after join; anchor your regular expressions. Each of these prevents a specific, common, silent error.
  • Compressed data should be streamed, not decompressed. BGZF plus tabix turns a compressed VCF into a randomly accessible one.
  • A pipeline’s exit status is the status of its last command only. set -euo pipefail at the top of every script converts silent failure into loud failure, which is the only kind you can fix.
  • A script that validates its arguments, writes errors to stderr and returns a meaningful exit status can be reused and trusted; one that does not is a liability inside a loop.
  • find with xargs or parallel scales a working command across many files; a workflow manager takes over when the steps acquire dependencies. On a cluster, measure what a job needs before requesting it, keep work off the login node, and use job arrays for per-sample work.

5.11 Exercises

  1. Using only grep, count how many sequences are in data/proteins/globins.fasta, and explain why grep -c '>' could give a different answer from grep -c '^>' on some FASTA files.
  2. Print the identifiers of the sequences in data/sequences/ that are longer than 10 kb, using bin/fasta_lengths.sh and one other command.
  3. How many records in data/annotation/mycoplasma-genitalium.gff are on the reverse strand? Use awk on column 7, restricting to gene features.
  4. Report the number of TP53 ClinVar records in each review-status category (CLNREVSTAT). Which is largest, and what does that imply for using the file clinically?
  5. Compute the mean read length and the GC content of data/reads/mgen_sim_R2.fastq.gz with a single awk program, and compare with the R1 figures given in this chapter.
  6. Write a script gc_report.sh that accepts one or more FASTA files and prints a three-column table of filename, total length and GC percentage. It must exit 2 with a usage message when given no arguments, and must skip unreadable files with a warning on stderr rather than stopping.
  7. The join in Section 5.5 produced 31 rows from 30 input genes. Modify the pipeline so that it reports one row per gene, and justify the rule you chose for resolving the duplicate.
  8. Compress data/variants/tp53_clinvar.vcf with gzip and with bgzip, compare the file sizes, then try to run tabix on each. Explain the result in terms of how the two formats store their blocks.
  9. A colleague’s script runs cat samples.txt | while read s; do process "$s" > "results/$s.txt"; done and reports success although three of 50 outputs are empty. Identify at least two defects and rewrite the loop.
  10. Write a SLURM submission script running exercise 6 as a job array over the FASTA files in data/. Explain what you would check with sacct afterwards.

5.12 Solutions

1. grep -c '^>' data/proteins/globins.fasta returns 4. Without the ^ anchor the pattern matches a > anywhere in the line. FASTA description lines frequently contain > inside a species name or a construct description, and some tools wrap descriptions onto continuation lines, so an unanchored count can exceed the true number of records.

2.

bin/fasta_lengths.sh data/sequences/mycoplasma-genitalium-NC_000908.2.fasta | awk '$2 > 10000 {print $1}'
NC_000908.2

Applying the same filter to every FASTA file at once gives both long sequences:

find data -name '*.fasta' -print0 | xargs -0 -n 1 bin/fasta_lengths.sh | awk '$2 > 10000 {print $1}'
NC_045512.2
NC_000908.2

3.

awk -F'\t' '$3 == "gene" {print $7}' data/annotation/mycoplasma-genitalium.gff | sort | uniq -c
    319 +
    227 -

227 of the 546 genes are on the reverse strand.

4.

bcftools query -f '%INFO/CLNREVSTAT\n' data/variants/tp53_clinvar.vcf | sort | uniq -c | sort -k1,1nr
   1731 criteria_provided,_single_submitter
   1145 criteria_provided,_multiple_submitters,_no_conflicts
    530 criteria_provided,_conflicting_classifications
    191 reviewed_by_expert_panel
    189 no_assertion_criteria_provided
     99 .
      3 no_classification_provided

The largest category is a single submitter providing assertion criteria: 1,731 of 3,888 records. Only 191 have been reviewed by an expert panel, 530 carry conflicting classifications and 99 have no value at all, printed by bcftools as .. Clinical significance in ClinVar is therefore not uniform evidence. Any use of the file for interpretation must filter on review status, the conflicting records must never be collapsed to a single label, and the missing-value category must be handled explicitly rather than silently dropped by a grep that happens not to match it.

5.

zcat data/reads/mgen_sim_R2.fastq.gz \
  | awk 'NR % 4 == 2 {total += length($0); gc += gsub(/[GCgc]/, ""); n++}
         END {printf "%d reads, mean length %.1f bp, GC %.2f%%\n", n, total / n, 100 * gc / total}'
20000 reads, mean length 150.0 bp, GC 31.74%

R2 contains the same 20,000 reads at the same 150 bp length as R1, and its GC content of 31.74% differs from R1’s 31.70% by less than the sampling noise you would expect from 3 Mb of sequence. Equal read counts in the two mates is the first thing to check in any paired library: an inequality means the files have been truncated or mismatched, and no aligner will tell you so politely.

6.

#!/usr/bin/env bash
# gc_report.sh - report length and GC content for one or more FASTA files.
set -euo pipefail

if [[ $# -eq 0 ]]; then
    echo "usage: $(basename "$0") FASTA [FASTA ...]" >&2
    exit 2
fi

for fasta in "$@"; do
    if [[ ! -r $fasta ]]; then
        echo "$(basename "$0"): cannot read '$fasta', skipping" >&2
        continue
    fi
    grep -v '^>' "$fasta" | tr -d '\n' \
      | awk -v name="$fasta" '{
            gc = gsub(/[GCgc]/, "")
            total = length($0) + gc
            printf "%s\t%d\t%.2f\n", name, total, 100 * gc / total
        }'
done
bin/gc_report.sh data/sequences/*.fasta
data/sequences/TP53-NM_000546.6.fasta   2512    53.38
data/sequences/mycoplasma-genitalium-NC_000908.2.fasta  580076  31.69
data/sequences/sars-cov-2-NC_045512.2.fasta 29903   37.97
bin/gc_report.sh data/sequences/sars-cov-2-NC_045512.2.fasta data/sequences/missing.fasta
bin/gc_report.sh
data/sequences/sars-cov-2-NC_045512.2.fasta 29903   37.97
gc_report.sh: cannot read 'data/sequences/missing.fasta', skipping
usage: gc_report.sh FASTA [FASTA ...]

Note that gsub returns the number of substitutions and destroys the G and C characters in $0, so the total length must be reconstructed as length($0) + gc.

7. Sorting by count and keeping the first occurrence of each key retains, for each gene symbol, the entry with the most GO terms. awk '!seen[$1]++' is the idiomatic first-occurrence filter and, unlike sort -u -k1,1, it is guaranteed to keep the row you sorted to the top:

join -t $'\t' work/network_genes.txt work/gene_go.tsv \
  | awk -F'\t' '{print $1 "\t" gsub(/GO:/, "")}' \
  | sort -k2,2nr | awk '!seen[$1]++' | wc -l
30

The justification has to be stated. The two CDKN2A entries carry 22 and 45 GO terms, and they exist because p16INK4a and p14ARF are translated from alternative reading frames of the same locus. If the analysis is gene-level — as a STRING network is — collapsing to the better-annotated entry is defensible provided it is documented. If the analysis is protein-level, the two entries are different proteins and must not be merged at all. Either way the decision belongs in the methods section, not in an undocumented sort -u.

8. gzip -c produces 88,847 bytes and bgzip -c produces 98,633 bytes from the same 1.1 MB VCF, so BGZF costs about 11% more space. Running tabix -p vcf on the plain gzip file fails:

[tabix] the compression of 'work/plain.vcf.gz' is not BGZF

Plain gzip is a single compressed stream: byte offset n of the compressed file cannot be located without decompressing everything before it, and restarting the decompressor mid-stream is impossible because the dictionary is cumulative. BGZF compresses the file in independent blocks of at most 64 kb and records a virtual offset for each, so the .tbi index can map a genomic coordinate to a block and the reader jumps straight to it. The 11% is what random access costs.

9. The defects are: the pipeline’s exit status reflects only the while loop, so a failing process is invisible; cat file | creates a subshell, so any variable set inside the loop is lost; read without -r mangles backslashes and strips leading whitespace; process failing after the shell has already created the output file leaves a zero-byte file that looks like a result; and nothing checks that the output is non-empty. A corrected version:

set -euo pipefail
while IFS= read -r sample; do
    [[ -n $sample ]] || continue
    if ! process "$sample" > "results/${sample}.txt.tmp"; then
        echo "failed: $sample" >&2
        rm -f "results/${sample}.txt.tmp"
        exit 1
    fi
    mv "results/${sample}.txt.tmp" "results/${sample}.txt"
done < samples.txt

Writing to a temporary name and renaming on success means a final output file exists only if the step that produced it succeeded — the same reasoning a workflow manager applies automatically.

10.

#!/usr/bin/env bash
#SBATCH --job-name=gc_report
#SBATCH --cpus-per-task=1
#SBATCH --mem=1G
#SBATCH --time=00:05:00
#SBATCH --array=1-6
#SBATCH --output=logs/gc_%A_%a.out
#SBATCH --error=logs/gc_%A_%a.err

set -euo pipefail

fasta=$(find data -name '*.fasta' | sort | sed -n "${SLURM_ARRAY_TASK_ID}p")
bin/gc_report.sh "$fasta"

Afterwards, sacct -j JOBID --format=JobID,JobName,Elapsed,MaxRSS,State shows elapsed time, peak memory and final state per task. MaxRSS is the number to act on: if it is a few megabytes, the 1 GB request was generous and can be reduced, which will shorten the queue wait next time. A task in state FAILED or OUT_OF_MEMORY points you at the matching logs/gc_*_*.err file.

5.13 Further reading

  • Wilson et al. (2017) is the shortest defensible statement of how to organise a computational project — file layout, naming, version control — and is worth reading before your next project rather than after it.
  • Sandve et al. (2013) gives ten rules for reproducible computational research; rules 1, 4 and 6 are directly about recording commands and intermediate results, which is what this chapter’s scripts are for.
  • Danecek et al. (2021) documents SAMtools and BCFtools, including bcftools query format strings and the BGZF/tabix indexing used in Section 5.6; it is the reference to keep open when working with VCF or BAM.
  • Danecek et al. (2011) is the original VCF specification paper, and is the place to confirm what each column and INFO key is allowed to contain before writing a parser of your own.
  • Köster and Rahmann (2012) and Di Tommaso et al. (2017) describe Snakemake and Nextflow, the two workflow managers most used in bioinformatics; read one of them when your shell loops acquire dependencies between steps.
  • Grüning et al. (2018) explains Bioconda, which is how most of the command-line tools in this book are installed reproducibly on a laptop or a cluster without administrator rights.
  • Ewels et al. (2016) describes MultiQC, the standard way to turn per-sample tool reports into one document; it pairs with the shell habits here because it consumes exactly the log files your pipelines produce.