Linux Command Line for Bioinformatics

Chapter 4 of 5: Pipes, redirection, loops and your first script

Pipes: join small tools into a pipeline

You have already used the pipe, |. It sends the output of one command straight into the next one, without a temporary file. This is the core idea of the command line: many small tools, each doing one job well, joined into a pipeline.

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

start zcat, awk and wc together
while zcat has data:
    zcat decompresses a chunk  -> writes it to the pipe
    awk reads lines from pipe  -> writes sequence lines to the next pipe
    wc reads lines             -> adds to its count
when zcat finishes, each program sees end-of-input and exits
wc prints the final count

Redirection: save results to files

grep -v '^#' data/variants.vcf > results/variants_body.txt
echo "analysis started" > results/log.txt
echo "variants extracted" >> results/log.txt
ls nosuchfile 2> results/errors.txt
cat results/log.txt results/errors.txt

• >: write output to a file; an existing file is replaced without warning
• >>: add output to the end of a file
• 2>: write error messages to a file
• <: read input from a file (we used it with tr in Chapter 2)
• > out.txt 2>&1: send both output and errors to the same file (order matters: file first, then 2>&1)

Exit codes: did the command work?

Every command finishes with an exit code: 0 means success, anything else means a problem. echo $? shows the code of the last command. You can chain commands on it: && runs the next command only if the first succeeded, and || runs it only if the first failed.

mkdir -p results/qc && echo "folder ready"
grep -q PASS data/variants.vcf && echo "has PASS variants" || echo "no PASS variants"

Loops: the same command for many files

Real projects have many samples. A for loop repeats commands for each item in a list:

for s in S01 S02 S03 S04; do
  echo "Processing $s"
done

The list can come from file names. The shell expands the pattern into matching names before the loop starts:

for f in data/genome.fa data/genes.gtf data/variants.vcf; do
  echo "$f: $(wc -l < "$f") lines"
done

In a real project the loop usually runs over samples, for example for r1 in raw/*_R1.fastq.gz; do … done, and calls a tool such as fastp or FastQC on each one. The FastQC tutorial uses exactly this pattern.

Your first script

When a set of commands works, save it in a script so you can run it again, share it and cite it. Create scripts/summary.sh with a text editor (nano scripts/summary.sh; save with Ctrl+O, Enter, exit with Ctrl+X) and type:

#!/usr/bin/env bash
# summary.sh: quick counts for the practice files
set -euo pipefail

genome="data/genome.fa"
reads="data/reads_R1.fastq.gz"
vcf="data/variants.vcf"

echo "Sequences in ${genome}: $(grep -c '^>' "$genome")"
echo "Reads in ${reads}: $(( $(zcat "$reads" | wc -l) / 4 ))"
echo "PASS variants in ${vcf}: $(awk -F'\t' '!/^#/ && $7 == "PASS"' "$vcf" | wc -l)"

Then make it executable and run it from the project folder:

chmod +x scripts/summary.sh
./scripts/summary.sh
./scripts/summary.sh > results/summary.txt

• #!/usr/bin/env bash: the first line, called the shebang, says which program runs the script
• set -euo pipefail: stop at the first error (-e), treat unset variables as errors (-u), and fail a pipeline if any part fails (pipefail)
• genome="data/genome.fa": a variable; no spaces around =, and use it as "$genome"
• chmod +x: give the file execute permission; ./ means “the script in this folder”

Find files and check disk space

find . -name '*.gtf'
du -sh data
df -h .

• find . -name '*.gtf': search this folder and everything below it by name
• du -sh data: how much space a folder uses
• df -h .: how much space is left on the disk you are on

Working on a server

Real datasets are analysed on a lab server, an HPC cluster or a cloud machine. You connect with ssh, copy files with rsync, and keep long jobs running with tmux or nohup. The server address and user name come from your lab or computing centre.

ssh [email protected]
rsync -avP data/ [email protected]:~/bhl-linux/data/
tmux new -s qc
nohup ./scripts/summary.sh > results/summary.log 2>&1 &

• ssh: log in to a remote machine; you get the same kind of terminal there
• rsync -avP: copy and resume interrupted transfers, keep file details (-a), show progress (-P); the trailing / on data/ copies the folder’s contents
• tmux new -s qc: start a session that keeps running if your connection drops; detach with Ctrl+B then D, come back with tmux attach -t qc
• nohup … &: run a command in the background and keep it running after you log out; its output goes to the log file