I'm working on trying to automate some BLAST searches. I need to pick up only the top three results from the BLAST results, however the parameter hitlist_size doesn't seem to be limiting my searches to only three results. No matter what size I specify I still end up getting > 3 hits for most samples (in others I get three, although I'm not sure if this is simply coincidence) . Does anyone have any insight to this?
import Bio
from Bio import SeqIO
from Bio.Blast import NCBIWWW
from Bio.Blast import NCBIXML
fasta_string = open("FASTA_Files/48_50.fasta").read()
result_handle = NCBIWWW.qblast("blastn", "nt", fasta_string,hitlist_size=3)
with open("XML_Files/48_50.xml", "w") as save_to:
save_to.write(result_handle.read())
result_handle.close()
The following FASTA sequence produces the expected number of hits for 46 but exceeds the expected number of hits for 47:
>46
NNNNNNNNNNGNNNNNNTGCAGTCGNANNNNNNNNNNNNNNNNNAGCTTGCTGCTTCGCTGACGAGTGGCGGACGGGTGA
GTAATGTCTGGGAAACTGCCTGATGGAGGGGGATAACTACTGGAAACGGTGGCTAATACCGCATAACGTCGCAAGACCAA
AGAGGGGGACCTTCGGGCCTCTTGCCATCAGATGTGCCCAGATGGGATTAGCTTGTTGGTGAGGTAACGGCTCACCAAGG
CGACGATCCCTAGCTGGTCTGAGAGGATGACCAGCCACACTGGAACTGAGACACGGTCCAGACTCCTACGGGAGGCAGCA
GTGGGGAATATTGCACAATGGGCGCAAGCCTGATGCAGCCATGCCGCGTGTATGAAGAAGGCCTTCGGGTTGTAAAGTAC
TTTCAGCGGGGAGGAAGGNGTTGTGGTTAATAACCGCAGCAATTGACGTTACCCGCANAANAAGCACCGGCTAACTCCGT
GCCAGCAGCCGCGGTAATACGGAGGGTGCANGCGTTAATCGNAATTACTGGNCGTAAAGCGCACGCAGGCGGTCTGTCAA
GTCTGATGTGAAATCCCCGGGCTCAACCTGGNAACTGCATTCNAAACTGGCAGGCTTGAGTCTTGNAGANGGGGGGTAGA
ATTCCAGGTGTANCGNTGAAATGCGTANAGATCTGGANGNAATACCGGTGGCNAANGCGGCCCCCTGGACAAAGACTGAC
GCTCANGTGCGAAAGCGTGGNGGAGCAAANAGGATTAGATACCCTGGTAGNCCNNCGCCNNANACGATGTCTACNTGGNA
GGNTTGTGNCCTTGAGGCGTGNCTNNCCNGNAGCTNAACGCGTTTAAGTANANCNNCCTGGGGCGAGCNACGGGCNGCNN
GGNTNAAAACNNNNNNTGNNNTTNNACGGGNGNNCCCCGCANNANCCGGCTGNNAGCATGNTGGATTNAANTTCGATNNN
NNCGCGAAGAANCCNNANNNNNGNNNNNNNANNNNNNNNNNAANNNTNNNNNANANNNNNNNNNNNNNNNGNNNNCTNNN
AGNANNGNNNCTGCATGGNNNNCNNNCNNGNTCNNGNNNNNNNNNNNNNGNNNNNNNNCCNNNANNNNNNNNNNNNTNAT
CNNNNNNNNNNNNNTNNNNNNNNNNNNNAGNNNNNNNNCNNNNNNNNNNANNTNNNAANN
>47
NNNNNNNNNNNNNNNNNNNNNNNNGNNNNNCGGGTNNCCNNNNGCTGGNTGNNNTGCTGACGAGTGGCGGACGGGTGAGT
AATGTCTGGGAAACTGCCTGAAGGGGGGGGATTCCTACTGGCCAGGGTGGCTAATACCGGGTAACGTCGNNNGANCAAAG
AGGGGGACCTTCGGGCCTCTTGCCATCACATGTGCCCGGATGGGATTAGCTTGTTGGTGAGGTAACGGCTCACCAAGGCG
ACGATCCCTAGCTGGTCTGAGAGGATGACCAGCCNCACTGGAACTGAGACACGGTCCCGACTCCTACGGGAGGCANCAGT
GGGGAATATTGCTCTTGGGCGCAAGCCTGATGCAGCCATGCCGCGGGTATGAGGAAGGCCTTCGGTTTGTAAAGTACTTT
CTCCGGGGAGGAAGGNGTNGTGGTGAATAACCGCTACANTTGANNCTNCCCGCNNAANAACCACCNGNTAACTCCNTGCN
NNNNGCCGCGGTAATACGGANGGTGCAAGNGTTAATCGNANTTACTGNNTGTTGAGCGCACGNNGGCGGCCTGTCNNNTC
TNATGTGAGATCCCCGGGCTCNCCCTGNNACCTGCATTCGNNNNNTGNNANGCTNGANTCTTGNNNNGNNGNGGNAGNAA
TTCCNNGTGTNNCGNNGAAATGCNNANAGATCTGNANANANNACNGGNGNCCAANGNNGNCCCCTGNTCTCNGACTGACG
CNNGAGTGCTGAANNGTGNAGAGCGNACAGGATTANANNNCCNGNTAGNCCGNCNCCNCACACCNATGTCTACATGNGAG
GNTNNNGNNNNNTGNGGCNNNNCNNTCCNNNAGCTNANGNNGTTNAANTANATCNNNCTNNNNCNNGCNNGGGGNCANNA
NGGGNNNAAANNTNNNNATNAATNTGACGGANNNNNCNNNNNNNNCNNNNNNANCATGNGGATTNANNNTNNNTNNNNNC
NNNNNANAACCNNANNNNNNNNNNNNNNTNNNNNNNANNNTNNNNNNNTNNNNNNGNNNNCNNNNNNNNACTNNNNNNNC
NNNNNNNNCANNGNNNNNNNNNNNANGNTNNNNNTGNNNNAANNNTGNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNTN
NNNNNNTNNNNNNNTNNNNTNNNNNNNNNNNNNNNNNNNNNNNTNNNNNNNNTCNNNNNN
The problem seems that a hit is something different than an alignment. A hit is a sequence match in the database, whereas an alignment is the actual position of the nucleotides. This could happen several times for the same sequence in the database.
In the case you tried, the xml that this saves has indeed three items under <Hit>, but the last hit you get has seven <Hsp> records, which seem to be what you are iterating through with the loop:
for rec in blast_records:
for alignment in rec.alignments:
for hsp in alignment.hsps:
print(rec.query) #DNA ID
Specifying alignments=3 instead of hits will get you 3 alignments maximum per hit, I think. You can specify both, if you want to control number of hits and number of alignments.