Summarize FASTQ and Annotation Files with Bash

Build one small Bash program that streams either FASTQ reads or a tab-separated annotation table, checks the limited assumptions described here, and writes a reproducible summary. The two input types need different checks, so you select the mode explicitly rather than guessing from the first line.

Prerequisites and tools

If you followed Part 2, continue in its already initialized shell environment; do not create a second Micromamba environment for this exercise. This summarizer uses Bash, awk, gzip, sort, and standard file utilities. It does not call SeqKit, R, or Python, so no additional packages are needed here. SeqKit is used for the sequence-aware steps in Part 2.

Before you run the examples

You should be comfortable with basic shell navigation and commands such as ls, cd, and cat. This script uses Bash-specific arrays and syntax, plus awk, gzip, sort, and standard file utilities. Run it with Bash, not sh. Check the required commands in your initialized environment:

bash --version | head -n 1
for tool in awk gzip sort cat mkdir dirname mktemp mv; do
  command -v "$tool" || exit 1
done
gzip --version | head -n 1

The FASTQ mode accepts the common four-line record layout: an identifier line beginning with @, a sequence, a separator beginning with +, and a quality string. It checks record structure and equal sequence/quality lengths. FASTQ variants and multiline records need a format-aware tool rather than this deliberately narrow summarizer. The annotation mode expects tab-separated columns in this order: chromosome, start, end, feature ID, strand. It calculates inclusive interval lengths as end - start + 1, so use one-based inclusive coordinates or adjust the formula for your file's coordinate convention.

Create a disposable example dataset

Use a fresh working directory. The setup stops rather than overwriting an existing shell-tutorial directory. It creates only synthetic teaching files, including one malformed FASTQ file for the safe failure check. Run this whole block from the directory where you want the lab folder created. The final cd leaves you inside it for the following steps.

[[ ! -e shell-tutorial ]] || { echo "shell-tutorial already exists; choose a fresh parent directory" >&2; exit 1; }
mkdir -p shell-tutorial/test_data shell-tutorial/results
cd shell-tutorial
: > .lesson-created

cat > test_data/annotations.tsv <<'TSV'
chr	start	end	feature_id	strand
chr1	100	200	geneA	+
chr1	300	350	geneB	-
chr2	50	400	geneC	+
TSV

cat > test_data/example.fastq <<'FASTQ'
@read1
ACTGACTGACTG
+
IIIIIIIIIIII
@read2
ACTGACTG
+
IIIIIIII
FASTQ
gzip -c test_data/example.fastq > test_data/example.fastq.gz

cat > test_data/bad-length.fastq <<'BADFASTQ'
@bad
ATGC
+
III
BADFASTQ

gzip -c writes compressed bytes to standard output and leaves the uncompressed fixture in place. The annotation table includes a header so the script can skip it explicitly. Do not use real, confidential, or irreplaceable data for this practice run.

Save the Bash summarizer

Here is the complete program. The quoted <<'BASH' delimiter makes Bash save the lines literally, without expanding variables or interpreting commands inside the program. Run this entire block from inside shell-tutorial. The chmod line makes the saved file executable; the later commands still invoke it with Bash explicitly.

The command form is ./summarize.sh {fastq|annotations} INPUT_FILE OUTPUT.tsv: choose the input layout, then give the input path and destination. Plain files and .gz files are supported, with compressed input integrity-tested before streaming. On success, an existing output at the same path is replaced, so choose a result filename you can safely overwrite.

[[ -f .lesson-created ]] || { echo "Run the setup first and stay inside shell-tutorial" >&2; exit 1; }
cat > summarize.sh <<'BASH'
#!/usr/bin/env bash
# Usage: summarize.sh {fastq|annotations} INPUT_FILE OUTPUT.tsv
set -euo pipefail

if [[ $# -ne 3 ]]; then
  echo "Usage: $0 {fastq|annotations} INPUT_FILE OUTPUT_TSV" >&2
  exit 2
fi
mode=$1
input=$2
output=$3
if [[ $mode != fastq && $mode != annotations ]]; then
  echo "Mode must be fastq or annotations" >&2
  exit 2
fi
if [[ $input != - && ! -r $input ]]; then
  echo "Input is not readable: $input" >&2
  exit 1
fi
if [[ $input != - && ( $input == "$output" || $input -ef $output ) ]]; then
  echo "Input and output must be different files: $input" >&2
  exit 1
fi
mkdir -p "$(dirname -- "$output")"
tmp=$(mktemp "${output}.tmp.XXXXXX")
trap 'rm -f -- "$tmp"' EXIT

if [[ $input == *.gz ]]; then
  gzip -t -- "$input"
  stream=(gzip -dc -- "$input")
elif [[ $input == - ]]; then
  stream=(cat)
else
  stream=(cat -- "$input")
fi

case $mode in
  fastq)
    {
      printf '#type\tcount\tmin_len\tmax_len\tmean_len\n'
      "${stream[@]}" | awk '
        { sub(/\r$/, "") }
        NR % 4 == 1 && substr($0, 1, 1) != "@" {
          print "Invalid FASTQ: record header must start with @ (line " NR ")" > "/dev/stderr"; bad=1; exit 1
        }
        NR % 4 == 2 {
          seq=$0
          if (seq == "") { print "Invalid FASTQ: empty sequence (line " NR ")" > "/dev/stderr"; bad=1; exit 1 }
          n=length(seq); sum+=n
          if (count==0 || n<min) min=n
          if (n>max) max=n
          count++
        }
        NR % 4 == 3 && substr($0, 1, 1) != "+" {
          print "Invalid FASTQ: separator line must start with + (line " NR ")" > "/dev/stderr"; bad=1; exit 1
        }
        NR % 4 == 0 && length($0) != length(seq) {
          print "Invalid FASTQ: sequence and quality lengths differ (record " count ")" > "/dev/stderr"; bad=1; exit 1
        }
        END {
          if (bad) exit 1
          if (NR == 0 || NR % 4 != 0) { print "Invalid FASTQ: empty or incomplete four-line record" > "/dev/stderr"; exit 1 }
          printf "reads\t%d\t%d\t%d\t%.2f\n", count, min, max, sum/count
        }'
    } > "$tmp"
    ;;
  annotations)
    {
      printf '#chr\tfeature_count\tavg_length\n'
      "${stream[@]}" | awk -F '\t' '
        NF==0 || $0 ~ /^#/ { next }
        tolower($1)=="chr" && tolower($2)=="start" { next }
        {
          if (NF < 3 || $1 == "" || $2 !~ /^[0-9]+$/ || $3 !~ /^[0-9]+$/ || $2 < 1 || $3 < $2) {
            print "Invalid annotation at line " NR ": expected chromosome and 1-based inclusive start/end with start <= end" > "/dev/stderr"; bad=1; exit 1
          }
          if (!($1 in count)) chromosomes++
          len=$3-$2+1; count[$1]++; sum[$1]+=len
        }
        END {
          if (bad) exit 1
          if (chromosomes==0) { print "No annotation rows found" > "/dev/stderr"; exit 1 }
          for (chr in count) printf "%s\t%d\t%.2f\n", chr, count[chr], sum[chr]/count[chr]
        }' | sort -k1,1
    } > "$tmp"
    ;;
esac
mv -f -- "$tmp" "$output"
trap - EXIT
echo "Wrote summary to $output" >&2
BASH
chmod 755 summarize.sh
ls -l summarize.sh

Run both modes and predict the results

Before running, predict the FASTQ summary: the two read lengths are 12 and 8, so the count is 2, the minimum is 8, the maximum is 12, and the mean is 10.00. For the annotation table, inclusive lengths are 101, 51, and 351 bases. That gives two chr1 features averaging 76.00 bases and one chr2 feature averaging 351.00 bases. These are arithmetic checks on invented inputs, not biological results.

Run the two commands from inside shell-tutorial:

./summarize.sh fastq test_data/example.fastq.gz results/fastq_summary.tsv
./summarize.sh annotations test_data/annotations.tsv results/annotations_summary.tsv
printf '\nFASTQ summary:\n'
cat results/fastq_summary.tsv
printf '\nAnnotation summary:\n'
cat results/annotations_summary.tsv

Expected table contents (the program also prints a short “Wrote summary” status message to standard error after each successful run):

#type	count	min_len	max_len	mean_len
reads	2	8	12	10.00

#chr	feature_count	avg_length
chr1	2	76.00
chr2	1	351.00

Try one safe failure path

The next command supplies the deliberately mismatched FASTQ. The program should reject it, remove its temporary output, and leave no summary file at the requested destination. This block treats the expected rejection as success; an unexpected successful run stops with an error.

if ./summarize.sh fastq test_data/bad-length.fastq results/rejected.tsv; then
  echo "Unexpectedly accepted malformed FASTQ" >&2
  exit 1
else
  echo "Rejected malformed FASTQ as expected"
fi
test ! -e results/rejected.tsv
if compgen -G 'results/rejected.tsv.tmp.*' > /dev/null; then
  echo "Temporary output was not cleaned up" >&2
  exit 1
fi

Expected diagnostic: Invalid FASTQ: sequence and quality lengths differ (record 1), followed by Rejected malformed FASTQ as expected. A malformed input must not be treated as a valid summary.

Modification exercise: add one short read

Predict the new summary before running the commands. Add a third six-base read, so the total length is 12 + 8 + 6 = 26 across three reads. The expected minimum is 6, maximum 12, and mean 8.67. The following appends a valid record to the invented uncompressed fixture, rebuilds its compressed copy, and writes a separate result so the first summary remains available:

cat >> test_data/example.fastq <<'FASTQ'
@read3
CCGTAA
+
IIIIII
FASTQ
gzip -c test_data/example.fastq > test_data/example-with-third.fastq.gz
./summarize.sh fastq test_data/example-with-third.fastq.gz results/fastq_summary_with_third.tsv
cat results/fastq_summary_with_third.tsv

Expected modified output:

#type	count	min_len	max_len	mean_len
reads	3	6	12	8.67

What the checks do and do not establish

• FASTQ mode rejects incomplete records, a misplaced header or separator, an empty sequence, or a sequence/quality length mismatch. A failed run removes its temporary output rather than installing the requested summary.
• A filename ending in .gz is integrity-tested with gzip -t and streamed through gzip. A failed integrity test stops the script. Compression format and suffix must agree; only a .gz suffix triggers decompression.
• Annotation rows with non-numeric coordinates, a start below 1, or an end before the start fail. Comment lines beginning with # and the described header row are skipped. The code does not validate chromosome names, strand values, feature IDs, overlap, genome build, or biological plausibility.
• Read length alone says nothing about per-base quality, adapter contamination, mapping, or whether trimming was appropriate. Use a dedicated sequencing QC workflow for those questions.

Why the shell checks matter

set -euo pipefail asks Bash to stop on many unhandled command failures, unset variables, and failed commands within a pipeline. In particular, pipefail prevents an upstream decompressor error from being hidden by a successful downstream command. The temporary file is created beside the destination and removed on failure; mv installs the completed file only after successful processing.

The awk program reads line records, checks the assumptions for its selected mode, and accumulates counts and lengths. In FASTQ mode, every fourth input line is treated as one record component, so wrapped or multiline sequence records are rejected rather than silently miscounted. Annotation output is sorted by chromosome name for stable row order. Bash arrays and [[ ... ]] are Bash features, so run the saved program with Bash rather than sh.

Common failures and safe adjustments

• gzip: ... not in gzip format or a failed integrity test: verify the file type and transfer integrity. Do not strip the suffix or process a damaged compressed file as plain text without checking it.
• FASTQ structure error: inspect the first records and confirm the file follows the four-line layout used here. FASTQ has variants; use a dedicated parser for wrapped records or other nonstandard layouts rather than weakening validation to accept ambiguous input.
• No annotation rows found: check that the table is tab-delimited, contains the stated columns, and has at least one data row after its header or comments. Change the explicit column mapping if your file uses a different schema.
• Lengths are off by one: confirm the source coordinate convention. The included end - start + 1 calculation is for 1-based inclusive intervals. For zero-based half-open intervals, length is end - start.
• Need to process many files: iterate over paths with a quoted loop rather than expanding a large wildcard into one command. From the directory containing the inputs, for example:

for file in data/*.fastq.gz; do
  [[ -e $file ]] || continue
  name=${file##*/}
  name=${name%.fastq.gz}
  ./summarize.sh fastq "$file" "results/${name}.summary.tsv"
done

The loop handles spaces in filenames and the no-match case; the script still checks each compressed input. Choose output names that do not collide with files you need to keep.

Permissions and longer-running jobs

chmod 755 summarize.sh gives the owner read, write, and execute permissions and gives group and other users read and execute permissions, subject to directory access controls. If others should not be able to read the script, choose a narrower mode such as chmod 700 summarize.sh.

For a long run, redirect both output streams to a log and run it in the background only when you understand the input size and destination:

nohup ./summarize.sh fastq large.fastq.gz results/large_summary.tsv > results/large_summary.log 2>&1 &

Record the job ID printed by the shell, inspect the log, and check completion before treating the output as finished. jobs lists jobs in the current shell; ps can inspect processes more broadly. A background process may still be running or may fail, so starting it is not evidence that the summary completed.

Optional cleanup

Cleanup is optional. Keep the lesson folder if you want to inspect the files later. To remove only this disposable lab, first return to the parent directory of shell-tutorial, confirm that it is the directory created above, and run the guarded command below. It refuses to remove a directory without the lesson marker. Do not put personal files in this disposable folder.

if [[ -f shell-tutorial/.lesson-created ]]; then
  rm -r -- shell-tutorial
else
  echo "Not removing shell-tutorial: lesson marker not found" >&2
  exit 1
fi

References

• Bash Reference Manual: pipelines, pipefail, redirection, and shell behavior.
• GNU Awk User’s Guide: records, fields, built-in functions, and output.
• GNU gzip manual: integrity testing and decompression to standard output.
• NCBI SRA file format guide: FASTQ structure and sequence/quality length relationship.
• UCSC data format FAQ: BED zero-based half-open coordinates and GFF one-based inclusive coordinates.

Scope and reproducibility

This is a small streaming summary, not a replacement for validated parsers, large-scale statistical analysis, or biological QC. It uses simple awk arithmetic and reports arithmetic mean lengths; for analyses needing other statistics, validated input semantics, or larger numeric ranges, use tools suited to that analysis. The local example is checked on the synthetic fixtures described above. Results can differ across awk and utility implementations, so test on the target system and record tool versions when portability matters. Do not run examples on irreplaceable or controlled data without following the relevant access and data-handling rules. Record input checksums and software versions when you need an auditable workflow.