AWK and sed on FASTA files

This page is optional. On it, you use AWK and sed on DNA sequences. You count sequences, see why a line in a sequence file is not the same as a sequence, and give the sequences shorter names without giving two of them the same name.

The sequences are 16S rRNA gene sequences of 11 bacterial strains. You use the same folder as on the previous pages, with the data in awk_data/.

Look at the files first

The sequences are stored in the FASTA format, one file per strain:

# The * stands for any text, so this lists all files ending in .fasta
ls awk_data/*.fasta

This prints:

awk_data/GCF_000009425.fasta
awk_data/GCF_000253315.fasta
awk_data/GCF_000714595.fasta
awk_data/GCF_014058685.fasta
awk_data/LjRoot303.fasta
awk_data/LjRoot44.fasta
awk_data/LjRoot59.fasta
awk_data/LjRoot60.fasta
awk_data/Root401.fasta
awk_data/Root68.fasta
awk_data/Root935.fasta

Look at the start of two of them:

head -3 awk_data/Root68.fasta
head -2 awk_data/GCF_000714595.fasta

This prints:

>Root68
AGAGTTTGATCATGGCTCAGATTGAACGCTGGCGGCAGGCCTAACACATGCAAGTCGAGCGGATGACAGGAGCTTGCTCC
TGAATTCAGCGGCGGACGGGTGAGTAATGCCTAGGAATCTGCCTGGTAGTGGGGGACAACGTTTCGAAAGGAACGCTAAT

>GCF_000714595-ECOLIN_RS01105_1
AAATTGAAGAGTTTGATCATGGCTCAGATTGAACGCTGGCGGCAGGCCTAACACATGCAAGTCGAACGGTAACAGGAATC

A FASTA file has two kinds of lines:

  • A header line starts with > and holds the name of the sequence.
  • The lines below it hold the sequence itself, until the next header. Here, each sequence is split over several lines of 80 letters.

The Root files hold one sequence each. The GCF files hold all 16S copies found in one genome, so they hold several sequences. Their names start with the genome ID, followed by a dash and more text.

Count the sequences

Each sequence has exactly one header, so counting the headers counts the sequences. In grep, as in sed, ^ means “the start of the line”, so '^>' finds lines that start with >. grep -c counts these lines:

# Per file
grep -c '^>' awk_data/*.fasta

# All files together
cat awk_data/*.fasta | grep -c '^>'

This prints:

awk_data/GCF_000009425.fasta:6
awk_data/GCF_000253315.fasta:6
awk_data/GCF_000714595.fasta:9
awk_data/GCF_014058685.fasta:7
awk_data/LjRoot303.fasta:1
awk_data/LjRoot44.fasta:1
awk_data/LjRoot59.fasta:1
awk_data/LjRoot60.fasta:1
awk_data/Root401.fasta:1
awk_data/Root68.fasta:1
awk_data/Root935.fasta:1

35

35 sequences in total, 9 of them from genome GCF_000714595.

ImportantAlways put quotes around ‘^>’

In the shell, > means “write the output into this file”. Without the quotes, grep -c ^> awk_data/Root68.fasta would not count anything. Instead the shell empties Root68.fasta to write into it, and grep waits for you to type something (stop it with Ctrl+C).

Line length is not sequence length

AWK can tell you how many characters a line has with length($0). But if you use it on a FASTA file, you get the length of each line, not of each sequence:

# The headers in this file: only one, so one sequence
grep '^>' awk_data/Root68.fasta

# For every line that does not start with >, print its length
awk '!/^>/ { print length($0) }' awk_data/Root68.fasta | head -3
awk '!/^>/ { print length($0) }' awk_data/Root68.fasta | tail -2
awk '!/^>/ { print length($0) }' awk_data/Root68.fasta | wc -l

This prints:

>Root68

80
80
80

80
26

18

Two new parts in this command:

  • /^>/ is a pattern that keeps lines that match the text between the slashes, here: lines that start with >. It works like grep '^>'.
  • ! means “not”. So !/^>/ keeps the lines that do not start with >, the sequence lines.

Root68 has one header, so it holds one sequence, but AWK prints 18 lengths: 17 lines of 80 letters and a last line of 26. The sequence is 17 × 80 + 26 = 1386 letters long.

length() is still useful on other files. It works on any text, also on a single column: length($1) gives the number of characters in column 1.

To get one length per sequence, you would have to add up the line lengths of each sequence. AWK can do that, but a tool made for sequence files does it in one command and leaves less room for mistakes. The tip below shows one.

seqkit is one example of a tool made for sequence files. seqkit fx2tab turns a FASTA file into a table. With -n it keeps only the names (not the sequences), and with -l it adds the length of each sequence:

# Name and length of one sequence
seqkit fx2tab -n -l awk_data/Root68.fasta

# Name and length of all sequences, saved as a table
seqkit fx2tab -n -l awk_data/*.fasta > lengths.tsv
wc -l lengths.tsv
head -3 lengths.tsv

This prints:

Root68  1386

35 lengths.tsv

GCF_000009425-LLMG_RS02645_19   1547
GCF_000009425-LLMG_RS10890_60   1547
GCF_000009425-LLMG_RS12185_65   1547

Root68 is 1386 letters long, the same as the sum of its line lengths above. To see the shortest and the longest sequence, sort by column 2. In sort -k2,2n, -k2,2 means “sort by column 2” and n means “as numbers”:

sort -k2,2n lengths.tsv | head -1
sort -k2,2n lengths.tsv | tail -1

This prints:

GCF_000714595-ECOLIN_RS27490_147    1363
GCF_000714595-ECOLIN_RS22205_97 1726

All 35 sequences are between 1363 and 1726 letters long. You need to install seqkit first, though. AWK is already on almost every system, and the same AWK patterns work on any text file, not only on sequences.

Shorter names with sed: spot the problem

Names like GCF_000714595-ECOLIN_RS01105_1 are long. For a figure, for example a tree, you may want shorter ones. The one thing you want to avoid when you change names: two sequences end up with the same name.

First, look at how the names are built. Each GCF name has two parts, separated by a dash: the genome ID (GCF_000714595), then a name for this one 16S copy within the genome (ECOLIN_RS01105_1). All 28 GCF names have exactly one dash. The 7 Root names have no dash, because each Root file holds only one sequence. So the dash is the place where you can cut a name, and each version below cuts in a different way.

Version 1: remove everything from the dash

In a regular expression, * means “the character before it, any number of times, or not at all”. So .* means “any text”, and -.* means “a dash and everything after it”:

cat awk_data/*.fasta | grep '^>' | sed 's/-.*//' | head -3

# Duplicate names
cat awk_data/*.fasta | grep '^>' | sed 's/-.*//' | sort | uniq -d

This prints:

>GCF_000009425
>GCF_000009425
>GCF_000009425

>GCF_000009425
>GCF_000253315
>GCF_000714595
>GCF_014058685

This version generates a problem: all 16S copies of one genome now have the same name, for example 9 times >GCF_000714595. The check shows this right away. This version is still useful if all you need is the genome ID, but not as new names for the sequences.

Version 2: from the dash to the underscore

This version keeps the number at the end: it replaces the dash and everything up to an underscore with a single underscore:

cat awk_data/*.fasta | grep '^>' | sed 's/-.*_/_/' | head -3

# Number of duplicate names
cat awk_data/*.fasta | grep '^>' | sed 's/-.*_/_/' | sort | uniq -d | wc -l

This prints:

>GCF_000009425_19
>GCF_000009425_60
>GCF_000009425_65

0

No duplicates, so the check passes. But look at the names: the whole middle part, LLMG_RS02645, is gone, not only LLMG. That is because .* takes as much text as it can. There are two underscores after the dash, and .*_ runs all the way to the last one. This is called a greedy match. A check that passes does not mean the result is what you want, so always look at the names too.

Version 3: from the dash to the first underscore

To stop at the first underscore, you need “any text, but without underscores”. Square brackets hold a list of characters, and a ^ right after the opening bracket turns the list around: [^_] means “any one character, except an underscore”. So -[^_]*_ means “a dash, then any characters that are not an underscore, then the first underscore”:

cat awk_data/*.fasta | grep '^>' | sed 's/-[^_]*_/_/' | head -3
cat awk_data/*.fasta | grep '^>' | sed 's/-[^_]*_/_/' | tail -3

# Number of duplicate names
cat awk_data/*.fasta | grep '^>' | sed 's/-[^_]*_/_/' | sort | uniq -d | wc -l

This prints:

>GCF_000009425_RS02645_19
>GCF_000009425_RS10890_60
>GCF_000009425_RS12185_65

>Root401
>Root68
>Root935

0

Only the letters between the dash and the first underscore (-LLMG) are gone, the rest of the name stays. The names stay unique, and the Root names, which have no dash, do not change.

QuestionExercise
  1. Use version 3 on all FASTA files together (cat awk_data/*.fasta), not only on the headers, and save the result as sequences_clean.fasta.
  2. Which checks would show that sed changed only the headers, and that the names are still unique? Think of three, then run them. Tip: cmp from the sed page compares two files.

Three things could go wrong: a line gets lost, two headers get the same name, or sed changes a sequence line. One check for each:

# 1. The sequence lines have no dash, so sed only changes the headers
cat awk_data/*.fasta | sed 's/-[^_]*_/_/' > sequences_clean.fasta

# 2a. The same number of lines
cat awk_data/*.fasta | wc -l
wc -l sequences_clean.fasta

# 2b. No duplicate names
grep '^>' sequences_clean.fasta | sort | uniq -d | wc -l

# 2c. Save the sequence lines before and after, and compare them
cat awk_data/*.fasta | grep -v '^>' > seq_before.txt
grep -v '^>' sequences_clean.fasta > seq_after.txt
cmp seq_before.txt seq_after.txt && echo 'identical'

This prints:

727
727 sequences_clean.fasta

0

identical

What you learned

  • A FASTA file has header lines starting with > and sequence lines below them. Count the sequences with grep -c '^>', always with the quotes.
  • /^>/ as an AWK pattern keeps the header lines, !/^>/ the other lines.
  • length($0) gives the length of a line, length($1) the length of column 1. A sequence may span several lines, so for sequence lengths use a tool made for sequence files, such as seqkit.
  • In a regular expression, .* is any text, and it takes as much as it can (greedy). [^_]* is any text without an underscore.
  • After shortening names, check for duplicates with sort | uniq -d, and look at the names yourself.

The next page, AWK on FASTQ files, is optional too. It uses AWK on sequencing reads. After it, Checking AWK or sed code from an LLM puts these checks together into a routine for any command you did not write yourself.