genetics/alg/blast.py

50 lines
1.8 KiB
Python
Raw Permalink Normal View History

2024-11-07 12:10:11 +00:00
from Bio.Blast import NCBIWWW, NCBIXML
from Bio import SeqIO
import os
# Ścieżka do pliku FASTA z sekwencją
fasta_file = "../data/sequence.fasta"
# Sprawdzenie, czy plik FASTA istnieje
if not os.path.exists(fasta_file):
print("Plik FASTA nie istnieje. Upewnij się, że plik sequence.fasta znajduje się w katalogu data.")
exit(1)
# Odczytywanie sekwencji z pliku FASTA
record = SeqIO.read(fasta_file, format="fasta")
sequence = str(record.seq)
# Wykonanie zapytania BLAST
print("Wysyłanie zapytania BLAST do NCBI...")
result_handle = NCBIWWW.qblast("blastn", "nt", sequence)
# Zapis wyników BLAST do pliku XML
output_file = "../data/blast_results.xml"
with open(output_file, "w") as out_handle:
out_handle.write(result_handle.read())
print(f"Wyniki BLAST zapisane do pliku {output_file}")
# Parsowanie wyników z pliku XML i wyświetlenie wyników
print("Analiza wyników BLAST...")
with open(output_file) as result_handle:
blast_records = NCBIXML.parse(result_handle)
for blast_record in blast_records:
for alignment in blast_record.alignments:
for hsp in alignment.hsps:
if hsp.expect < 0.01: # wartość e-value
print("\n****HIT****")
print(f"Sequence: {alignment.title}")
print(f"Length: {alignment.length}")
print(f"E-value: {hsp.expect}")
print(f"Score: {hsp.score}")
print(f"Identities: {hsp.identities}/{hsp.align_length}")
print("Query sequence:")
print(hsp.query[0:75] + "...")
print("Match:")
print(hsp.match[0:75] + "...")
print("Subject sequence:")
print(hsp.sbjct[0:75] + "...")