Parsing FASTA with SeqIO

A FASTA file is the most common way to store DNA, RNA, or protein sequences as plain text, where each entry has a header line starting with a greater-than sign followed by the sequence letters. You could split the file by hand, but real files wrap sequences across many lines and can contain blank lines, so manual parsing easily corrupts your data. Biopython is a free Python library for biology, and its SeqIO reader parses FASTA files correctly for you.

>seq1 Homo sapiens beta-globin coding start fragment
ATGGTGCACCTGACTCCTGAGGAGAAGTCTGCCGTTACT
>seq2 Homo sapiens insulin exon fragment
ATGGCCCTGTGGATGCGCCTCCTGCCCCTGCTGGCGCTG
>seq3 Escherichia coli lacZ start fragment
ATGACCATGATTACGGATTCACTGGCCGTCGTTTTACAACGTCGTGACTGGGAAAACC

Save these six lines in a plain text file named sequences.fasta in the same folder as your Python script. On each header line the greater-than sign marks the start, the first word such as seq1 is the identifier, and the rest of the line is a free-text description. The letters A, T, G, and C on the lines below are the DNA bases that make up that sequence.

Reading with SeqIO

from Bio import SeqIO

for record in SeqIO.parse("sequences.fasta", "fasta"):
    print("ID:", record.id)
    print("Description:", record.description)
    print("Length:", len(record.seq))
    print("Sequence:", record.seq)
    print()

ID: seq1
Description: seq1 Homo sapiens beta-globin coding start fragment
Length: 39
Sequence: ATGGTGCACCTGACTCCTGAGGAGAAGTCTGCCGTTACT

ID: seq2
Description: seq2 Homo sapiens insulin exon fragment
Length: 39
Sequence: ATGGCCCTGTGGATGCGCCTCCTGCCCCTGCTGGCGCTG

ID: seq3
Description: seq3 Escherichia coli lacZ start fragment
Length: 58
Sequence: ATGACCATGATTACGGATTCACTGGCCGTCGTTTTACAACGTCGTGACTGGGAAAACC

The first line imports the SeqIO reader from Biopython. SeqIO.parse takes the file name and the format name written as fasta, then hands you one SeqRecord object per entry as the loop turns, so a huge file never has to fit in memory all at once. Each record gives you record.id (the first word of the header), record.description (the full header text after the greater-than sign), and record.seq, a Seq object whose length you measure with the built-in len function.

Counting and summarizing

from Bio import SeqIO

# list() collects every SeqRecord so we can count and reuse them
records = list(SeqIO.parse("sequences.fasta", "fasta"))

# build one length number per record
lengths = [len(record.seq) for record in records]

print("Number of records:", len(records))
print("Lengths:", lengths)
print("Shortest:", min(lengths))
print("Longest:", max(lengths))
print("Total bases:", sum(lengths))
print("Average length:", sum(lengths) / len(records))

Number of records: 3
Lengths: [39, 39, 58]
Shortest: 39
Longest: 58
Total bases: 136
Average length: 45.333333333333336