Linux Command Line for Bioinformatics

Chapter 3 of 5: Search and count in FASTA, FASTQ, GTF and VCF files

Most everyday questions in bioinformatics are simple counts: how many sequences, how many reads, how many genes, how many variants passed the filter? This chapter answers each one with a short command. Run them inside ~/bhl-linux.

FASTA: count sequences and measure their lengths

In a FASTA file, every sequence starts with a header line beginning with >, followed by one or more sequence lines. So the number of sequences equals the number of lines that start with >.

grep -c '^>' data/genome.fa
grep '^>' data/genome.fa

• grep pattern file: print lines that match the pattern
• -c: print only the number of matching lines
• '^>': ^ means “start of the line”, so only lines that begin with > match

Always put the pattern in quotes. Without quotes, the shell reads > as “write output to a file”. The command grep -c > data/genome.fa would not search anything; it would empty your genome file.

To measure each sequence, we use awk, a small language for line-by-line processing. This one-liner remembers the name at each header and adds up the lengths of the sequence lines below it:

awk '/^>/ {if (name) print name, len; name = substr($1, 2); len = 0; next}
     {len += length($0)}
     END {print name, len}' data/genome.fa

for each line in the file:
    NR = line number; $0 = the line; $1, $2 ... = its fields
    if the line starts with ">":            # pattern /^>/
        if a name is stored: print name, len
        name = first field without ">"      # substr($1, 2)
        len = 0
        skip to the next line               # next
    otherwise:
        len = len + length of the line
after the last line:                        # END
    print name, len

GC content

GC content is the share of G and C bases among all bases: GC% = 100 × (G + C) ÷ (A + C + G + T). Two ways to get it:

grep -v '^>' data/genome.fa | tr -d '\n' | fold -w1 | sort | uniq -c

awk '!/^>/ {gc += gsub(/[GCgc]/, ""); at += gsub(/[ATat]/, "")}
     END {printf "GC = %.1f%%\n", 100*gc/(gc+at)}' data/genome.fa

Check it by hand: (79 + 64) ÷ (81 + 79 + 64 + 61) = 143 ÷ 285 = 50.2%. The first command is easy to read but slow on big genomes (one line per base); the awk version is the one to keep.

FASTQ: count reads the right way

A FASTQ read always takes exactly four lines: header, sequence, a + line, and quality. So the number of reads is the number of lines divided by four.

echo $(( $(zcat data/reads_R1.fastq.gz | wc -l) / 4 ))

• $( ... ): run the command inside and use its output as text
• $(( ... )): do whole-number arithmetic

The line number also tells you what each line is. NR % 4 is the remainder after dividing the line number by 4: it is 2 for every sequence line and 0 for every quality line.

zcat data/reads_R1.fastq.gz | awk 'NR % 4 == 2' | head -n 2
zcat data/reads_R1.fastq.gz | awk 'NR % 4 == 2 {print length($0)}' | sort -n | uniq -c

The second command is a read-length histogram: all 6 reads are 36 bases long. After adapter trimming (see our FastQC and MultiQC tutorial) you would see several lengths here.

grep can also search the reads. The Illumina adapter starts with AGATCGGAAGAGC, and an N marks a base the machine could not call:

zcat data/reads_R1.fastq.gz | grep -c AGATCGGAAGAGC
zcat data/reads_R1.fastq.gz | awk 'NR % 4 == 2' | grep -c N

GTF and VCF: work with columns

GTF (genes) and VCF (variants) are tab-separated tables. Each line is one record, and each column is a field you can pick by number.

Count features in a GTF file

cut -f3 data/genes.gtf | sort | uniq -c

• cut -f3: keep column 3 (the feature type); cut splits on tabs by default
• sort: put identical lines next to each other
• uniq -c: collapse neighbouring identical lines and count them

previous = none; count = 0
for each line:
    if line == previous: count = count + 1
    else:
        if previous is not none: print count, previous
        previous = line; count = 1
print count, previous

awk can filter rows and do arithmetic on columns. Here we list each gene with its chromosome, start, end, strand and length:

awk -F'\t' '$3 == "gene" {print $1, $4, $5, $7, $5 - $4 + 1}' data/genes.gtf

-F'\t' tells awk that columns are separated by tabs. This matters: by default awk splits on any run of spaces and tabs, so the GTF attribute column, which contains spaces, would be cut into extra fields.

Count and filter variants in a VCF file

A VCF file starts with header lines beginning with #. grep -v inverts the match, so it keeps only the variant lines:

grep -v '^#' data/variants.vcf | wc -l
grep -v '^#' data/variants.vcf | cut -f7 | sort | uniq -c

Now the two best PASS variants by quality (column 6), sorted from high to low:

awk -F'\t' '!/^#/ && $7 == "PASS"' data/variants.vcf | sort -t$'\t' -k6,6nr | cut -f1,2,4,5,6 | head -n 2

• !/^#/ && $7 == "PASS": not a header line and FILTER equals PASS
• sort -t$'\t' -k6,6nr: sort by column 6 only, numerically (n), in reverse (r)
• cut -f1,2,4,5,6: keep CHROM, POS, REF, ALT and QUAL

Which command for which file

Defaults that trip people up