Counting & filtering sequences

In earlier lessons you met FASTA files, the plain-text format biologists use to store DNA and protein sequences. Each sequence, also called a record, begins with a header line that starts with a greater-than sign (>), followed by the letters of the sequence itself. In this lesson you will learn to count those records and pull out just their headers, all by typing short commands into the shell (the text window where you type instructions to your computer).

Finding lines with grep

grep is a small program that reads through a file and prints only the lines that match a pattern you give it. The pattern '^>' means 'lines that begin with >', because ^ stands for the start of a line and > is the character we want; the single quotes keep the shell from meddling with those special symbols. Later we will add a flag to make grep count matches instead of printing them.

grep '^>' reads.fasta

>read1 Homo sapiens
>read2 Homo sapiens
>read3 Escherichia coli
>read4 Escherichia coli

Counting sequences

grep -c '^>' reads.fasta

4

Because every record has exactly one header line, counting the headers is the same as counting the sequences, and the -c flag reports that count as a single number. A pipe, written as the vertical bar (|), takes the output of one command and feeds it into the next, so command1 | command2 lets them work as a team. awk is a tiny language for working with columns: it splits each line into fields wherever it sees spaces or tabs, so $1 means the first field (here, the read name right after the >).

grep '^>' reads.fasta | awk '{print $1}'

>read1
>read2
>read3
>read4