genetics/alg/blast.py

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 = "../out/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] + "...")