50 lines
1.8 KiB
Python
50 lines
1.8 KiB
Python
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] + "...")
|
|
|