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" >sequence_1 ATGCTGCTGCTCGTCGTC >sequence_2 TTTATATATAAGTGTAGTAGA >sequence_3 GGGCGCGCGAGATAATATATTAT
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)
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") print(tmp_sequence) >>> ATGCTGCTGCTCGTCGTC # fetch part of the sequence "sequence 1" from position 0 to 3 tmp_sequence = sequences_object.fetch("sequence_1", 0, 3) print(tmp_sequence) >>> 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!