Command-line text processing for bioinformatics

Use line-oriented shell commands for quick inspection and simple tables. For sequence records, first check whether the format and operation can safely be handled one line at a time.

Learning objectives and prerequisites

Prerequisites: Basic shell navigation, a Bash terminal, the usual text utilities, and SeqKit. Start in a fresh practice directory, not a directory containing research data: the setup writes the five named teaching files. The table exercises need only the shell utilities; complete the SeqKit preflight before the sequence-aware exercises.

Learning objectives: After working through this lesson, you will be able to:

• Use grep (global regular expression print), sed (stream editor), cut, sort, uniq, and pipes to filter and reshape biological text files.
• Write short, reproducible one-line or small multi-step shell workflows that extract specific records and produce simple summaries.
• Recognize when line-oriented tools are appropriate versus when sequence-aware tools are required.

Setup checkpoint before the exercise

Use a Bash terminal, not an R or Python console. On Windows, use a Linux/WSL terminal for this Bioconda route rather than native PowerShell. You do not need R, Python, Scanpy or Seurat for this lesson. Check the tools before creating the practice files.

Already have SeqKit? Run this readiness check. It stops at the first missing command without closing your interactive terminal. If every command is found and a SeqKit version prints, skip the installation steps and continue to the exercise.

(
  for tool in awk cut sort uniq grep sed seqkit; do
    command -v "$tool" || { printf 'Missing tool: %s\n' "$tool" >&2; exit 1; }
  done
  seqkit version
)

Need SeqKit? First check whether Micromamba is available:

command -v micromamba && micromamba --version

If it is not found, stop here and follow only Install Micromamba and verify the executable, then return to this checkpoint. You do not need to complete that later lesson's QC environment lab or its optional R/Python examples first. On a managed computer, use your institution's approved installation route.

Once the executable is available, load its Bash hook in this terminal. This enables activation for the current shell; it does not rewrite your startup files. Stop if either command fails:

hook="$(micromamba shell hook --shell bash)" && eval "$hook"

For the first setup only, create a separate environment containing SeqKit. Review the package transaction and confirm it when prompted. Do not install it into a research environment. If cli-sequences already exists from an earlier visit, skip creation and try the activation check instead; do not overwrite it.

micromamba create -n cli-sequences --override-channels -c conda-forge -c bioconda --channel-priority strict seqkit

After creation succeeds, activate it and verify the actual command. The chained checks do not print a success from another installed SeqKit if activation fails:

micromamba activate cli-sequences &&
command -v seqkit &&
seqkit version

Rerun the readiness check above before continuing. Expected evidence is a command path for each utility and a SeqKit version, not just a downloaded file or an environment name. In a new terminal, repeat the hook and activation checks, not environment creation. If solving fails, retain the error and inspect platform support and channel order instead of switching channels blindly. See the official Micromamba installation guidance and Bioconda channel policy.

Setup validation: On Ubuntu Linux x86_64, the displayed check, Bash hook, creation and activation commands succeeded with Micromamba 2.9.0 and SeqKit 2.14.0 in a fresh home directory and environment root, reusing an existing executable and package cache. The readiness check and both expected FASTA motif matches passed. The installer itself and other operating systems were not retested in this revision.

Why text processing matters in bioinformatics

Many bioinformatics data formats (FASTA, FASTQ, variant tables, gene annotation TSVs) are plain text. Being able to filter, extract, and summarize these files from the command line lets you explore data quickly, create reproducible preprocessing steps, and plug into larger pipelines. These tools process input a line at a time and can be useful in pipelines over large files; sorting may need substantial memory or temporary disk space, so do not assume every command in a pipeline has constant memory use. That same line-oriented behavior means you must know the file layout (for example, whether sequence lines are wrapped) before running a command like grep or cut.

Key concept: line orientation vs. record orientation

Tools such as grep, sed, cut, sort, and uniq work on lines. If a biological record spans multiple lines (common in wrapped FASTA), naive line-based commands can yield incorrect results. Sequence-aware tools (for example, seqkit) treat sequences as records and should be used when format complexity matters.

Reproducible workflow: filter, extract, and summarize

Below is a sequential workflow demonstrating typical tasks. Each code block lists the expected input files and what each command produces. Use exact paths if your files are in another directory (for example, /home/user/data/sequences.fasta).

Files used in this example

• sequences.fasta: three invented records, including one whose motif is split across two sequence lines.
• reads.fastq: two invented eight-base reads, with equal sequence and quality lengths. malformed.fastq deliberately has four bases but only three quality characters.
• data.tsv: a tab-delimited table with a header line and multiple columns (e.g., gene_id, sample, count).

Run the setup in a fresh practice area. These sequences, IDs, counts and quality characters are invented solely for this lesson; they are not biological observations. To rerun setup while inside cli-part2-practice, first run cd .. to return to its parent directory and check pwd; otherwise the setup creates a nested practice directory. Re-running from the same parent overwrites the five teaching filenames. For comparison or recovery, the same inputs are available as FASTA, valid FASTQ, malformed FASTQ, TSV, and ID list.

mkdir -p cli-part2-practice
cd cli-part2-practice
printf 'gene_id\tsample\tcount\nGENE_A\tS1\t120\nGENE_B\tS1\t45\nGENE_A\tS2\t90\nGENE_C\tS2\t200\n' > data.tsv
printf 'GENE_A\nGENE_C\n' > genes_of_interest.txt
cat > sequences.fasta <<'FASTA'
>gene_A_unwrapped forward_motif
ATGCGTAA
>gene_B_wrapped split_motif
CCA
TGCGTT
>gene_C_negative no_forward_motif
TTTTAAAA
FASTA
cat > reads.fastq <<'FASTQ'
@read_1
ATGCGTAA
+
IIIIIIII
@read_2
CCATGCGT
+
IIIIIIII
FASTQ
cat > malformed.fastq <<'BAD_FASTQ'
@bad_length
ATGC
+
III
BAD_FASTQ
head -n 2 data.tsv

The first two TSV lines should be its header and GENE_A/S1/120, separated by tabs. Inspect the fixtures before processing them. The split sequence CCA followed by TGCGTT is one record, not two reads.

1) Inspect files safely

Before running extraction commands, inspect the file to determine layout.

# Show the first 40 lines of a FASTA to check header and wrap behavior
head -n 40 sequences.fasta

# View the file interactively (use q to quit)
less sequences.fasta

What these produce: text printed to your terminal showing the file start. Look for whether sequence headers begin with ">" (FASTA) and whether sequences continue on multiple lines.

2) Count sequence records in a FASTA file (safe line-based and record-aware)

If each FASTA header begins with ">" and headers are on separate lines, counting headers is a quick proxy for record count. If sequences are wrapped, this remains valid because headers are still single lines.

# Count FASTA headers (sequences.fasta -> counts to stdout)
grep -c '^>' sequences.fasta

Predict, then check: the supplied FASTA should give 3. This counts header lines, not sequence lines. It is a record count only under the usual FASTA convention of one header per record; it does not validate the sequences, their provenance, or ID uniqueness.

If you prefer a sequence-aware count using seqkit (recommended for complex FASTA files):

# Requires seqkit installed; replace path or install in your environment
seqkit stats sequences.fasta

SeqKit should report three FASTA records and a total length of 25 bases for this fixture. Inspect errors and the input format; successful parsing cannot establish that the file represents the intended biological data.

3) Extract headers and tidy them (remove the leading >)

# Extract header lines and remove the leading '>'
grep '^>' sequences.fasta | sed 's/^>//' > headers.txt
# Output: headers.txt contains one header per line (header text only)

Explanation: grep selects header lines; sed substitutes the leading character. headers.txt is ready to be used for lookups or counts.

4) Compare a line hit with a complete record

Predict: will a line search and a sequence-record search find the same two motif-containing records? Read sequences.fasta before running either command.

# Line search: includes the following line, not the preceding record header
# The next line may actually be the NEXT record's header.
grep -n -A 1 'ATGCG' sequences.fasta > matching_lines.txt

# Record search: join wrapped sequence lines and search only the forward strand
seqkit grep -s -P -p 'ATGCG' sequences.fasta > matching_records.fasta
seqkit seq -n -i matching_records.fasta

Worked result: the line search finds only ATGCGTAA on line 2; its extra context line is the next record's header, not the matching record's header. SeqKit returns gene_A_unwrapped and gene_B_wrapped. The second motif crosses CCA and TGCGTT, so it cannot appear in a single-line match.

-s searches sequence content; -P restricts this exercise to the forward strand. Without -P, SeqKit also searches the reverse complement by default. Neither a motif match nor a matching record count establishes function or experimental validity.

5) Convert the toy FASTQ, then challenge the assumption

The supplied valid FASTQ has exactly four lines per record. The following AWK example is a demonstration for that fixture only, not a validator or a general FASTQ parser. It copies the header and sequence while discarding quality information.

# Demonstration ONLY for the supplied valid four-line teaching fixture
awk 'NR%4==1{printf ">%s\n", substr($0,2)} NR%4==2{print}' reads.fastq > reads.toy.fasta

Expect two FASTA records, read_1 and read_2, each eight bases long. A plausible output does not prove the original FASTQ was valid: the AWK command never checked its plus lines or sequence/quality length agreement.

For general FASTQ conversion, use the record-aware parser. It can fail after reading part of a file, so this block writes to a temporary file and publishes the final filename only after a successful exit. This is format handling, not biological quality control.

# Run in a child Bash shell so a parse failure stops before publishing output.
bash <<'BASH'
set -euo pipefail
input=reads.fastq
output=reads.checked.fasta
tmp=$(mktemp "${output}.tmp.XXXXXX")
trap 'rm -f "$tmp"' EXIT
seqkit fq2fa "$input" > "$tmp"
mv "$tmp" "$output"
BASH

For the valid fixture, reads.checked.fasta should match reads.toy.fasta. Quality characters are intentionally lost in both outputs; keep the original FASTQ when downstream steps need quality scores.

Deliberate failure: malformed.fastq has four bases and three quality characters. Predict what each method will do, then run:

# This creates plausible-looking FASTA even though the supplied FASTQ is broken.
# It is a counterexample, NOT a conversion method for unknown inputs.
awk 'NR%4==1{printf ">%s\n", substr($0,2)} NR%4==2{print}' malformed.fastq > misleading.fasta

# The record-aware parser must reject this malformed fixture.
if seqkit fq2fa malformed.fastq > rejected.fasta.tmp; then
  rm -f rejected.fasta.tmp
  printf 'Unexpected acceptance: stop and inspect the parser/input.\n' >&2
  false
else
  rm -f rejected.fasta.tmp
  printf 'Rejected the malformed FASTQ as expected.\n'
fi

Worked result: AWK writes >bad_length and ATGC and exits successfully. SeqKit rejects the unequal sequence and quality lengths with a nonzero exit; the exact numeric error code can vary. The expected final message is Rejected the malformed FASTQ as expected. Do not use partial output from a failed parse, and do not treat parser acceptance as proof of correct samples, quality encoding or biological provenance.

6) Work with tab-delimited data (column selection, filtering, counting)

# Show header and first 5 rows (safe inspection)
head -n 6 data.tsv

# Extract columns 1 and 3 to a new file (tab-delimited)
cut -f 1,3 data.tsv > col1_col3.tsv
# Output: col1_col3.tsv contains only columns 1 and 3

# Filter rows where column 3 (count) > 100 (assumes header present)
# This keeps the header and filters numeric rows
(head -n1 data.tsv && awk -F'\t' 'NR>1 && $3>100' data.tsv) > filtered_counts.tsv
# Output: filtered_counts.tsv with header plus rows where col3 > 100

# Count column-1 values (excluding the header); the first line must be a header
cut -f1 data.tsv | tail -n +2 | sort | uniq -c | sort -k1,1nr > gene_counts.txt
# Output: column-1 value frequencies in descending count order (ties are not explicitly ordered)

Notes: The parentheses around the header-preserve pipeline ensure the header is copied, then rows are filtered. Replace column indices based on your file schema. With the supplied table, the count filter produces GENE_A/S1/120 and GENE_C/S2/200. Column-1 frequencies are 2 for GENE_A and 1 each for GENE_B and GENE_C (the tied values may appear in either order).

7) Small example: extract records for a set of gene IDs

The setup block created genes_of_interest.txt with two IDs that occur in the sample table. First-column matching matters: a whole-line search can match an ID in another field or as part of a longer string.

# genes_of_interest.txt contains gene IDs, one per line
# Fixed-string substring search across each entire line; not a column-exact match
grep -F -f genes_of_interest.txt data.tsv > genes_subset_grep.tsv
# Exact match against column 1; FILENAME test also works for an empty ID list
awk -F'\t' 'FILENAME == ARGV[1] {a[$1]=1; next} ($1 in a)' genes_of_interest.txt data.tsv > genes_subset.tsv
# Output: three matching data rows; empty if there are no matches

Inspect with cat genes_subset.tsv; the expected first fields are GENE_A, GENE_A, and GENE_C. To test no matches, temporarily replace the list with NOT_PRESENT; the output should contain zero rows. Restore the original list afterward. This result intentionally omits the TSV header.

Output inspection and interpretation

What to look for after each step:

• Counts from grep -c '^>' or seqkit stats: confirm the number matches expectations (e.g., sample metadata). A mismatch can indicate extra header-like lines, concatenated files, or file corruption.
• Headers.txt contents: verify that IDs are what downstream tools expect (no extra whitespace, no > character).
• Matching_records.fasta: check that full records were preserved; if only partial sequences appear, you may have used a line-based grep on wrapped FASTA.
• Filtered TSV outputs: confirm the header presence and that numeric comparisons (such as $3 > 100) behave as intended. Do not rely on an AWK numeric comparison to validate a field: nonnumeric text can be coerced or compared in ways that do not express the intended numeric threshold. Validate the field and normalize it before filtering.

What a plausible result means and what it does NOT prove:

• A list of sequence IDs matching 'ATGCG' indicates those sequences contain that motif. It does NOT imply functional significance or experimental validity; further validation and context are required.
• A high read count or many matching records may indicate biological signal or contamination; text processing alone cannot distinguish these possibilities.
• Counting headers gives the number of records in the file, but not the uniqueness of IDs. Use uniq or additional checks to test for duplicated IDs.

Try a variation and explain the result

1. Change the forward-strand motif from ATGCG to TTTT. Predict the matching ID before running the record-aware command.
2. Repeat the checked-conversion block in a fresh practice directory with input=malformed.fastq and output=must-not-publish.fasta. Predict the exit status and whether that final file should exist.

What was checked

The displayed fixture setup, table filters and ID matching, wrapped-motif comparison, both valid conversions and malformed-record rejection were run on invented inputs under Ubuntu Linux x86_64 with Bash 5.2.21, mawk 1.3.4 (20240123), GNU grep 3.11, GNU Coreutils 9.4, and SeqKit 2.14.0. The failed checked-conversion variation was tested for temporary-file cleanup and no new final output. The Micromamba preparation route, other operating systems, performance on large files and biological analyses were not tested here.

Troubleshooting

Below are common symptoms, likely causes, and suggested fixes.

Symptom: grep returns no matches, but you expect hits

• Cause: pattern case mismatch. Fix: use -i for case-insensitive search: grep -i 'pattern' file.
• Cause: pattern contains characters interpreted by the shell. Fix: quote patterns properly (single quotes preferred) and escape regex metacharacters, or use grep -F for fixed-string matching.
• Cause: line wrapping or multi-line records. Fix: use a record-aware tool like seqkit or reformat the file so each record is single-line if appropriate.

Symptom: cut results include empty columns or wrong fields

• Cause: different delimiter than tab. Fix: identify delimiter (use head | sed or less) and supply -d option to cut, or use awk -F to specify delimiter.
• Cause: header contains extra whitespace. Fix: clean headers with sed to trim whitespace: sed 's/^[ \t]*//; s/[ \t]*$//' file.

Symptom: awk numeric comparison fails

• Cause: non-numeric characters (commas, units) in the numeric column. Fix: remove or convert non-numeric characters before numeric comparison, e.g., gsub(/,/, "", $3) inside awk.

Symptom: pipelines run slowly on very large files

• Cause: using algorithms that are not streaming-friendly or too many intermediate writes to disk. Fix: keep operations piped where possible (avoid unnecessary temporary files), or use specialized tools (seqkit, samtools for SAM/BAM) and appropriate hardware.

Limitations, version and reference-data caveats, and responsible use

• Line-based tools assume a consistent file structure. If input formats vary across files, automated scripts can silently produce incorrect outputs. Always inspect representative samples with head or less.
• Command behavior can vary slightly with tool versions (GNU vs BSD variants). When writing reproducible pipelines, document tool versions (for example, grep --version) and consider using a package manager (conda/mamba) to pin versions.
• On very large datasets, file encoding and newline styles (Windows CRLF vs Unix LF) can affect pattern matching. Normalize files with dos2unix if needed.
• Responsible use: extracting sequence IDs or motifs does not replace proper ethical review or consent checks for sensitive genomic data. Do not share sensitive outputs without appropriate permissions.
• Record-aware tools (seqkit, samtools, bcftools) are preferable for formats with multi-line records or binary representations (e.g., BAM/CRAM). Prefer these when available rather than attempting complex sed/awk workarounds.

Next-step checklist

• Verify file layout for each dataset (wrapped vs unwrapped FASTA, tab delimiter).
• Create and test one small, documented command-line workflow for a common task (e.g., extract gene set, count reads).
• Record tool versions used and consider containerizing workflows (Docker/Singularity) or using an environment manager (conda) for reproducibility.
• When uncertain about format complexities, use seqkit or other format-aware tools rather than relying solely on line-based utilities.
• Practice safe data handling and check institutional policies before sharing processed data.

References

• GNU grep manual: documentation and options for pattern searching.
• GNU sed manual: stream editor reference for substitution and simple text transformations.
• GNU Coreutils (cut, sort, uniq, wc): authoritative reference for standard text-processing utilities.
• SeqKit usage documentation: a fast toolkit for FASTA/Q manipulation with record-aware operations.
• FASTA format guidance (NCBI): definition-line identifiers and sequence-format requirements.