AWK on FASTQ files

This page is optional. On it, you use AWK on sequencing reads. You count the reads in a file, see why the obvious grep command counts too many, learn to pick out one line of each read with AWK, and turn the reads into a FASTA file.

You use the same folder as on the previous pages, with the data in awk_data/.

Look at the files first

The reads come from an amplicon sequencing run, stored in the FASTQ format. There are two files: a sample that worked well, and a sample that failed and has only a few reads:

ls awk_data/*.fastq.gz

This prints:

awk_data/Run1_S137_R1.fastq.gz
awk_data/Run1_S158_R1.fastq.gz

The .gz at the end means that the files are compressed, as sequencing files almost always are. To read them, use zcat: it unpacks the file and prints the text, but leaves the file itself as it is. You know zcat from the CLI workshop. On a Mac, use gzcat instead, here and in the rest of this page.

Look at the first read of the small file:

zcat awk_data/Run1_S158_R1.fastq.gz | head -4

This prints:

@M08325:104:000000000-LT57L:1:2115:28775:11445 1:N:0:TGATGTAAGT
GACCCTCAGCAATCTTAAACTTCTTAGACGAATCACCAGAACGGAAAACATCCTTCATAGAAATTTCACGCGGCGGCAAGTTGCCATACAAAACAGGGTCGCCAGCAATATCGGTATAAGTCAAAGCACCTTTAGCGTTAAGGTACTGAATCTCTTTAGTCGCAGTAGGCGGAAAACGAACAAGCGCAACAGTAAACATAGTGCCATGCTCAGGAACAAAGAAACGCGGCACAGAATGTTTATAGGTCTGTTGAACACGACCAGAAAACTGGCCTAACGACGTTTGGTCAGTTCCATGAAA
+
C<C<<EGGFDCFGGGGGGFGGEFFGGGFB,88AF<?4@,EDACEDFGCFFCFFEAFFFCCFFF9CFEGGFFCFCEGFGGCFGGGFGGDFFGGGGGFGEFGGCFGGGGGGGGGGGGGGGFFFF9FGFGGFDF<DFFFFCCFD9<FFGC,C=@FFDDFFE9ECCEE7CFC,>>*@EGGGCCEGGGC88EGG*8,<9AFG7E?:EDC7EA*?+;C475=FCFF@CGFG>>57BD<*)<?GCFF:65<>BB44(5@*@D?B:>2),:.4575:F<)7733@8)394461(-.4A<4<:<9<)--)

Each read always takes 4 lines: a header that starts with @, the sequence, a line with +, and a quality line with one character per base. Each character stands for how sure the sequencer is about that base. The CLI workshop explains the four lines in more detail.

Count the reads

Because every read takes exactly 4 lines, you can count the lines and divide by 4:

zcat awk_data/Run1_S158_R1.fastq.gz | wc -l

This prints:

40 lines, so 10 reads.

On the FASTA page, you counted sequences with grep -c '^>'. Every header in a FASTQ file starts with @, so grep -c '^@' looks like the same idea. Before you run it, predict: what number should it print?

zcat awk_data/Run1_S158_R1.fastq.gz | grep -c '^@'

This prints:

One read too many, and no error. To see where the extra one comes from, let grep -n print the line number in front of each line it finds. cut -c 1-40 keeps only the first 40 characters of each line, so the output stays readable:

zcat awk_data/Run1_S158_R1.fastq.gz | grep -n '^@' | cut -c 1-40

This prints:

1:@M08325:104:000000000-LT57L:1:2115:287
5:@M08325:104:000000000-LT57L:1:2115:217
9:@M08325:104:000000000-LT57L:1:2106:121
13:@M08325:104:000000000-LT57L:1:2102:24
17:@M08325:104:000000000-LT57L:1:1119:19
21:@M08325:104:000000000-LT57L:1:1115:47
25:@M08325:104:000000000-LT57L:1:1113:27
29:@M08325:104:000000000-LT57L:1:1113:55
32:@CBCCGCFGDGFCF:FGGGGGGC8B@8BFFGGGCGGG
33:@M08325:104:000000000-LT57L:1:1109:22
37:@M08325:104:000000000-LT57L:1:1105:22

The headers are on lines 1, 5, 9 and so on, every 4 lines. Line 32 breaks the pattern. It is not a header, but the quality line of read 8. @ is one of the characters that stand for a quality, so a quality line can start with @, and grep cannot tell it apart from a header.

In a file with 10 reads, this is one read in ten. In a large file, the error is smaller, but it is there in almost every file, and it changes from file to file. You see this in the exercise below.

Pick one line of each read: `NR % 4`

You can count the reads with AWK instead, by looking at the position of a line, not at how it starts. The header is always line 1 of a read: line 1, 5, 9, and so on. AWK finds these lines with NR % 4 == 1.

% gives the remainder after dividing: NR % 4 is what is left of NR after you take away as many 4s as you can. For line 5, that is 1 (5 = 4 + 1). For line 8, it is 0 (8 = 4 + 4). Print it next to the first 8 lines to see the pattern. cut -c 1-30 again keeps the lines short:

zcat awk_data/Run1_S158_R1.fastq.gz | head -8 | cut -c 1-30 | awk '{ print NR, NR % 4, $0 }'

This prints the line number, NR % 4, and the start of the line:

1 1 @M08325:104:000000000-LT57L:1:
2 2 GACCCTCAGCAATCTTAAACTTCTTAGACG
3 3 +
4 0 C<C<<EGGFDCFGGGGGGFGGEFFGGGFB,
5 1 @M08325:104:000000000-LT57L:1:
6 2 AATATCTATAACAACTATTTTAAAGCGCCG
7 3 +
8 0 CCCCCGFGG9FFGGGGGGGGGGGGGG7FGF

The remainder repeats 1, 2, 3, 0 for each read. So:

Line of the read	`NR % 4`
header	1
sequence	2
`+`	3
quality	0 (not 4)

Now count the headers:

zcat awk_data/Run1_S158_R1.fastq.gz | awk 'NR % 4 == 1' | wc -l

This prints:

The same as 40 lines divided by 4. You can also ask AWK directly which quality lines start with @. /^@/ is a pattern for lines that start with @, as /^>/ was on the FASTA page:

zcat awk_data/Run1_S158_R1.fastq.gz | awk 'NR % 4 == 0 && /^@/ { print NR }'

This prints:

The same line that grep counted as a header.

NR % 4 relies on one thing: every read takes exactly 4 lines. That is true for the FASTQ files you get from current sequencers, but check it once on a new file: the number of lines should divide by 4 without a remainder.

AWK on a .gz file reads junk, without an error

If you forget zcat, AWK reads the compressed bytes as if they were text:

awk 'NR % 4 == 1' awk_data/Run1_S158_R1.fastq.gz | wc -l

This prints:

One “header” instead of 10, and no error. Without wc -l, AWK prints unreadable characters to your screen. For any .gz file, start with zcat file | ....

From FASTQ to FASTA

Some tools want FASTA, not FASTQ. A FASTA file needs only the header, with > instead of @, and the sequence. So keep line 1 and line 2 of each read:

zcat awk_data/Run1_S158_R1.fastq.gz | awk 'NR % 4 == 1 { print ">" substr($0, 2) } NR % 4 == 2 { print }' > Run1_S158_R1.fasta

# Check: 2 lines per read, and 10 headers
wc -l Run1_S158_R1.fasta
grep -c '>' Run1_S158_R1.fasta
head -2 Run1_S158_R1.fasta | cut -c 1-60

This prints:

20 Run1_S158_R1.fasta
10
>M08325:104:000000000-LT57L:1:2115:28775:11445 1:N:0:TGATGTA
GACCCTCAGCAATCTTAAACTTCTTAGACGAATCACCAGAACGGAAAACATCCTTCATAG

Taken apart:

NR % 4 == 1 { print ">" substr($0, 2) }: on the header line, print a > and then the header without its first character. substr($0, 2) is the part of the line from character 2 to the end, so the @ is dropped.
NR % 4 == 2 { print }: on the sequence line, print the line as it is.
The + line and the quality line match neither pattern, so AWK leaves them out.

Here grep -c '>' is a safe check, because > never appears in a sequence line.

Exercise

Count the reads in Run1_S137_R1.fastq.gz, the sample that worked well, in two ways: with grep -c '^@', and with NR % 4 == 1. Before you run them, predict: will the two numbers be the same?
Which check tells you which of the two numbers is right?

Click me to see an answer

zcat awk_data/Run1_S137_R1.fastq.gz | grep -c '^@'
zcat awk_data/Run1_S137_R1.fastq.gz | awk 'NR % 4 == 1' | wc -l

# Check: the number of lines, divided by 4
zcat awk_data/Run1_S137_R1.fastq.gz | wc -l

This prints:

2489
2485
9940

9940 lines divided by 4 is 2485 reads, the same as the AWK count. grep counts 4 reads too many. You can find the 4 quality lines that start with @:

zcat awk_data/Run1_S137_R1.fastq.gz | awk 'NR % 4 == 0 && /^@/ { print NR }'

This prints:

Here the error is less than 0.2%, easy to miss. In the small file it was 10%. You find it only if you compare with a second count.

You may also see grep -c '^@M08325' in code from others, with the start of the header instead of only @. In this file it gives 2485 as well, because all reads come from the same sequencer, M08325. In another file, the header starts with something else, so look at your headers first, and still compare with the number of lines divided by 4.

What you learned

A FASTQ file holds 4 lines per read: header, sequence, +, quality. Read a .fastq.gz file with zcat. AWK on the compressed file gives junk, without an error.
grep -c '^@' can count too many reads, because a quality line can also start with @. Count the lines and divide by 4, or use NR % 4 == 1.
% gives the remainder after dividing. NR % 4 == 1 picks the header of each read, == 2 the sequence, == 0 the quality line.
substr($0, 2) is the line from character 2 on. With it, NR % 4 == 1 { print ">" substr($0, 2) } NR % 4 == 2 { print } turns FASTQ into FASTA.

On the next page, Checking AWK or sed code from an LLM, you put these checks together into a routine for any command you did not write yourself.

Look at the files first

Count the reads

Pick one line of each read: NR % 4

From FASTQ to FASTA

What you learned

Pick one line of each read: `NR % 4`