BLAST with Biopython

from Bio.Blast import NCBIWWW

# Use a known accession query from Biopython BLAST examples
# (this tends to return many hits more reliably than a very short random sequence)
query = "8332116"

# Submit an online BLASTN search to NCBI and ask for XML output
result_handle = NCBIWWW.qblast(
    program="blastn",
    database="nt",
    sequence=query,
    hitlist_size=5,
    format_type="XML",
)

# Save results for later parsing steps
with open("blast_result.xml", "w", encoding="utf-8") as out_handle:
    out_handle.write(result_handle.read())

result_handle.close()
print("Saved BLAST results to blast_result.xml for query:", query)

from Bio.Blast import NCBIXML

# Parse the saved BLAST XML file
with open("blast_result.xml", "r", encoding="utf-8") as result_file:
    blast_record = NCBIXML.read(result_file)

print("Query:", blast_record.query)
print("Alignments returned:", len(blast_record.alignments))

# Show summary for the first three hits
for alignment in blast_record.alignments[:3]:
    hsp = alignment.hsps[0]
    print("Hit title:", alignment.title)
    print("E-value:", hsp.expect)
    print("Alignment length:", hsp.align_length)
    print("-" * 60)

from Bio.Blast import NCBIXML

# Load parsed BLAST results
with open("blast_result.xml", "r", encoding="utf-8") as result_file:
    blast_record = NCBIXML.read(result_file)

# Keep only strong matches by e-value threshold
threshold = 1e-10
similar_hits = []

for alignment in blast_record.alignments:
    for hsp in alignment.hsps:
        if hsp.expect <= threshold:
            similar_hits.append(
                {
                    "title": alignment.title,
                    "evalue": hsp.expect,
                    "identity_fraction": hsp.identities / hsp.align_length,
                }
            )
            break

print("Strong similar hits found:", len(similar_hits))
for hit in similar_hits[:5]:
    print(hit)

from Bio.Blast import NCBIXML

# Parse BLAST XML
with open("blast_result.xml", "r", encoding="utf-8") as result_file:
    blast_record = NCBIXML.read(result_file)

# Biopython BLAST records can expose query length as query_letters
query_len = getattr(blast_record, "query_length", None) or getattr(blast_record, "query_letters", None)
if not query_len:
    raise ValueError("Could not determine query length from BLAST record.")

ranked = []
for alignment in blast_record.alignments:
    best_hsp = max(alignment.hsps, key=lambda h: h.identities / h.align_length)
    query_coverage = best_hsp.align_length / query_len
    identity_fraction = best_hsp.identities / best_hsp.align_length
    ranked.append(
        {
            "title": alignment.title,
            "evalue": best_hsp.expect,
            "identity_fraction": identity_fraction,
            "query_coverage": query_coverage,
        }
    )

# Sort by identity first, then coverage
ranked.sort(key=lambda x: (x["identity_fraction"], x["query_coverage"]), reverse=True)

if not ranked:
    print("No BLAST hits were returned for ranking.")
else:
    for hit in ranked[:5]:
        print(hit)

import csv
from Bio.Blast import NCBIXML

# Parse BLAST XML
with open("blast_result.xml", "r", encoding="utf-8") as result_file:
    blast_record = NCBIXML.read(result_file)

query_len = getattr(blast_record, "query_length", None) or getattr(blast_record, "query_letters", None)
if not query_len:
    raise ValueError("Could not determine query length from BLAST record.")

rows = []
for alignment in blast_record.alignments:
    hsp = alignment.hsps[0]
    rows.append(
        {
            "title": alignment.title,
            "evalue": hsp.expect,
            "identities": hsp.identities,
            "align_length": hsp.align_length,
            "identity_fraction": hsp.identities / hsp.align_length,
            "query_coverage": hsp.align_length / query_len,
        }
    )

fieldnames = [
    "title",
    "evalue",
    "identities",
    "align_length",
    "identity_fraction",
    "query_coverage",
]

with open("blast_hits_summary.csv", "w", newline="", encoding="utf-8") as csvfile:
    writer = csv.DictWriter(csvfile, fieldnames=fieldnames)
    writer.writeheader()
    if rows:
        writer.writerows(rows)

print(f"Wrote blast_hits_summary.csv with {len(rows)} hit rows")

import time
from Bio.Blast import NCBIWWW

queries = {
    "query_1": "8332116",
    "query_2": "AY851612",
}

for query_id, sequence in queries.items():
    for attempt in range(1, 4):
        try:
            handle = NCBIWWW.qblast(
                program="blastn",
                database="nt",
                sequence=sequence,
                hitlist_size=3,
                format_type="XML",
            )
            with open(f"{query_id}.xml", "w", encoding="utf-8") as f:
                f.write(handle.read())
            handle.close()
            print(f"Saved {query_id}.xml")
            break
        except Exception as exc:
            print(f"{query_id} attempt {attempt} failed: {exc}")
            if attempt == 3:
                raise
            time.sleep(5)
    # Respectful delay between NCBI requests
    time.sleep(2)