Just today I wrote a script to do exactly this using Biopython
. I also add a warning If any the headers I wanted to filter was not present in the original fasta. So the python script filter_fasta_by_list_of_headers.py
is :
#!/usr/bin/env python3
from Bio import SeqIO
import sys
ffile = SeqIO.parse(sys.argv[1], "fasta")
header_set = set(line.strip() for line in open(sys.argv[2]))
for seq_record in ffile:
try:
header_set.remove(seq_record.name)
except KeyError:
print(seq_record.format("fasta"))
continue
if len(header_set) != 0:
print(len(header_set),'of the headers from list were not identified in the input fasta file.', file=sys.stderr)
and usage :
filter_fasta_by_list_of_headers.py input.fasta list_of_scf_to_filter > filtered.fasta
P.S. it's quite easy to turn over the script to extract the sequences from the list (just the print line would have to move after line header_set.remove(seq_record.name
)
I made a modification to use in
and benchmarked it on normal insect genome assembly - 1Gbp, 300,000 scaffolds and I gave it lists of 5,000 and 100,000 scaffolds to filter. Regardless of the size of the list or algorithm used it took 29-30s. In other words, if LBYL is faster than EAFP, it's on a scale that totally does not matter.