The Fastest Way to Read a FASTA in Python



This tutorial teaches a fast approach to how to read sequences from large FASTA files in Python using Pysam. I know how reading large FASTA files can be painful, so I hope this tutorial is helpful for you.

1. Reading a large FASTA File Very Quickly

When working with a FASTA file in Python, I load all the sequences into a dictionary where the sequences id is the key, the sequence is the value, and I use the dictionary to run my analysis.

However, in many cases, loading all the sequences into a dictionary is not the most efficient approach to work with large FASTA files, for example, when fetching sequences from the NR/NT database.

Now, Pysam can create a FASTA object and fetch the whole sequence given the sequence ID or a slice of it.

Please assume the following FASTA file with three sequences:

# file name "test.fasta"

Using the subsequent pysam calls, it creates a FASTA object for your input file

from pysam import FastaFile

fasta = "test.fasta"
# read FASTA file
sequences_object = FastaFile(fasta)

When calling “FastaFile,” Pysam calls “samtools faidx” which indexes your FASTA file if not present.

2. Fetching Regions of a Sequence from a FASTA file

Once the FASTA is indexed, it guarantees the agile FASTA reading and fetching. We can read the whole sequence or part of it as demonstrated below.

# fetch whole sequence
tmp_sequence = sequences_object.fetch("sequence_1")


# fetch part of the sequence "sequence 1" from position 0 to 3
tmp_sequence = sequences_object.fetch("sequence_1", 0, 3)

>>> ATG

Last but not least, here, we used a small FASTA file as a test. Nonetheless, when using a large FASTA file (500 Gbp) is, when Pysam fetching method highlights the most.

I hope that this tutorial was helpful for you. Reading large FASTA files can be painful, and I’m sure the approach presented here will help you read it quickly and make the task less painful; enjoy it!

Related Posts