#! /usr/bin/env python3 # vim:fenc=utf-8 # # Copyright © 2024 user # # Distributed under terms of the MIT license. import os from Bio import Entrez, SeqIO from Bio.Seq import Seq from Bio.SeqFeature import SeqFeature, FeatureLocation from Bio.Restriction import RestrictionBatch, EcoRI, BamHI, HindIII, Bpu1102I from Bio.Restriction import AllEnzymes, Analysis from Bio.Graphics import GenomeDiagram from reportlab.lib import colors from Bio.Blast import NCBIWWW, NCBIXML # Ustawienia Entrez Entrez.email = "baiobelfer@gmail.com" # Wpisz swój adres email # Funkcja do pobierania sekwencji z NCBI i zapisywania jej lokalnie def fetch_and_save_sequence(db, id, file_path): handle = Entrez.efetch(db=db, id=id, rettype="gb", retmode="text") record = SeqIO.read(handle, "genbank") handle.close() SeqIO.write(record, file_path, "genbank") return record # Funkcja do wczytywania sekwencji z lokalnego pliku def load_local_sequence(file_path): with open(file_path, 'r') as file: record = SeqIO.read(file, "genbank") return record #%% # ID dla sekwencji w NCBI gen_id= "4e46" plasmid_id = "pET-28a(+)" # Ścieżki do lokalnych plików gen_file_path = '../data/fasta/'+gen_id+'.fasta' plasmid_file_path = '../data/gb/'+plasmid_id+'.gb' #%% # Sprawdzanie i pobieranie sekwencji genu if not os.path.exists(gen_file_path): print(f"Plik {gen_file_path} nie istnieje. Pobieranie z serwera NCBI...") gen_record = fetch_and_save_sequence("nuccore", gen_id, gen_file_path) print(f"Pobrano i zapisano sekwencję genu do pliku {gen_file_path}.") else: print(f"Plik {gen_file_path} istnieje. Wczytywanie danych z lokalnego pliku...") gen_record = load_local_sequence(gen_file_path) print(f"Wczytano dane z pliku {gen_file_path}.") gen_seq = str(gen_record.seq) print(f"Gen Sequence: {gen_seq[:60]}...") #%% # Sprawdzanie i pobieranie sekwencji plazmidu if not os.path.exists(plasmid_file_path): print(f"Plik {plasmid_file_path} nie istnieje. Pobieranie z serwera NCBI...") record = fetch_and_save_sequence("nuccore", plasmid_id, plasmid_file_path) print(f"Pobrano i zapisano sekwencję plazmidu do pliku {plasmid_file_path}.") else: print(f"Plik {plasmid_file_path} istnieje. Wczytywanie danych z lokalnego pliku...") record = load_local_sequence(plasmid_file_path) print(f"Wczytano dane z pliku {plasmid_file_path}.") #%% # Funkcja do wczytania sekwencji z pliku FASTA def load_fasta_sequence(file_path): with open(file_path, "r") as fasta_file: record = SeqIO.read(fasta_file, "fasta") return record #%% # Wczytanie sekwencji z pliku FASTA fasta_record = load_fasta_sequence(gen_file_path) print(f"Wczytano sekwencję: {fasta_record.id}") #%% # Krok 2: Wykonanie BLAST do identyfikacji RefSeq def blast_and_find_refseq(sequence_record): print("Wykonywanie BLAST...") result_handle = NCBIWWW.qblast("blastp", "refseq_protein", sequence_record.seq) # Zapisz wynik BLAST w formacie XML with open("blast_result.xml", "w") as out_handle: out_handle.write(result_handle.read()) result_handle.close() # Parsowanie wyników BLAST with open("blast_result.xml") as result_handle: blast_record = NCBIXML.read(result_handle) # Pobranie najlepszego dopasowania RefSeq best_hit = blast_record.alignments[0] refseq_id = best_hit.accession print(f"Znaleziono RefSeq ID: {refseq_id}") return refseq_id # Krok 3: Pobranie pliku GenBank dla zidentyfikowanego RefSeq ID def fetch_and_save_genbank(db, id, file_path): print(f"Pobieranie sekwencji dla ID: {id} z bazy danych {db}") with Entrez.efetch(db=db, id=id, rettype="gb", retmode="text") as handle: record = SeqIO.read(handle, "genbank") SeqIO.write(record, file_path, "genbank") print(f"Sekwencja została zapisana w pliku {file_path}") # Wyszukiwanie RefSeq dla sekwencji z pliku FASTA refseq_id = blast_and_find_refseq(fasta_record) # Ścieżka do zapisu pliku GenBank genbank_file_path = '../data/gb/' + refseq_id + '.gb' # Pobranie i zapisanie pliku GenBank na podstawie RefSeq ID if refseq_id: fetch_and_save_genbank("protein", refseq_id, genbank_file_path) else: print("Nie znaleziono odpowiedniego RefSeq ID.") #%% # Wyświetlanie etykiet (annotacji) dla CDS for feature in record.features: if feature.type == "CDS": print(f"Type: {feature.type}") print(f"Location: {feature.location}") # Pobierz numer pozycji nukleotydów dla CDS start = feature.location.start end = feature.location.end strand = '+' if feature.location.strand == 1 else '-' print(f"Start: {start}, End: {end}, Strand: {strand}") # Pobierz nazwę genu, jeśli jest dostępna if 'gene' in feature.qualifiers: print(f"Gene: {feature.qualifiers['gene'][0]}") # Pobierz nazwę produktu białkowego, jeśli jest dostępna if 'product' in feature.qualifiers: print(f"Product: {feature.qualifiers['product'][0]}") # Wyświetl translację, jeśli jest dostępna if 'translation' in feature.qualifiers: print(f"Translation: {feature.qualifiers['translation'][0]}") print("\n" + "-"*50 + "\n") #%% # Importujemy bibliotekę BioPython from Bio import SeqIO # Przechodzimy przez wszystkie funkcje w rekordzie i wypisujemy informacje for feature in record.features: if feature.type in ["CDS", "promoter", "terminator", "misc_feature", "RBS", "protein_bind"]: print(f"Type: {feature.type}") print(f"Location: {feature.location}") print(f"Start: {feature.location.start}") print(f"End: {feature.location.end}") print(f"Strand: {'+' if feature.strand == 1 else '-'}") # Kierunek nici if 'label' in feature.qualifiers: print(f"Label: {feature.qualifiers['label'][0]}") if 'gene' in feature.qualifiers: print(f"Gene: {feature.qualifiers['gene'][0]}") if 'product' in feature.qualifiers: print(f"Product: {feature.qualifiers['product'][0]}") if 'note' in feature.qualifiers: print(f"Note: {feature.qualifiers['note'][0]}") print() # Dla czytelności #%% # Wyświetl etykiety (annotacje) for feature in record.features: print(f"Type: {feature.type}") print(f"Location: {feature.location}") if 'gene' in feature.qualifiers: print(f"Gene: {feature.qualifiers['gene']}") if 'product' in feature.qualifiers: print(f"Product: {feature.qualifiers['product']}") print() #%% # Przejdź przez wszystkie funkcje w rekordzie i wypisz ich szczegółowe informacje for feature in record.features: if feature.type in ["promoter", "RBS","CDS", "protein_bind", "misc_feature", "rep_origin", "terminator"]: print(f"Type: {feature.type}") print(f"Location: {feature.location}") print(f"Strand: {'+' if feature.strand == 1 else '-'}") # Szczegółowe opisy if feature.type == "CDS": gene_name = feature.qualifiers.get('gene', ['Unknown gene'])[0] product = feature.qualifiers.get('product', ['Unknown product'])[0] print(f"Gene Name: {gene_name}") print(f"Product: {product}") elif feature.type == "protein_bind": binding_site = feature.qualifiers.get('bound_moiety', ['Unknown binding site'])[0] print(f"Binding Site: {binding_site}") elif feature.type == "misc_feature": note = feature.qualifiers.get('note', ['No additional information'])[0] print(f"Note: {note}") elif feature.type == "rep_origin": print("This is a replication origin.") # Pobierz sekwencję funkcji feature_seq = record.seq[feature.location.start:feature.location.end] print(f"Sequence ( ): {feature_seq}\n") # Uwzględnij orientację nici (strand) if feature.strand == -1: feature_seq = feature_seq.reverse_complement() print(f"Sequence (^): {feature_seq}\n") # Znajdź i wyświetl sekwencję terminatora T7 for feature in record.features: if feature.type == "terminator" and "T7" in feature.qualifiers.get('note', [''])[0]: print("Terminator T7 znaleziony:") print(f"Type: {feature.type}") print(f"Location: {feature.location}") print(f"Sequence: {record.seq[feature.location.start:feature.location.end]}") print() #%% # Znajdź i wyświetl sekwencję promotora T7 sq = [] for feature in record.features: if feature.type == "promoter" and "T7" in feature.qualifiers.get('note', [''])[0]: print("Promotor T7 znaleziony:") print(f"Type: {feature.type}") print(f"Location: {feature.location}") # Pobierz sekwencję promotora promoter_seq = record.seq[feature.location.start:feature.location.end] # Sprawdź orientację i odwróć/uzupełnij jeśli potrzebne if feature.location.strand == 1: strand_orientation = "(+)" sequence_with_ends = f"5'-{promoter_seq}-3'" elif feature.location.strand == -1: strand_orientation = "(-)" promoter_seq = promoter_seq.reverse_complement() sequence_with_ends = f"3'-{promoter_seq}-5'" else: strand_orientation = "(unknown strand)" sequence_with_ends = str(promoter_seq) #print(f"Strand: {strand_orientation}") print(f"Sequence: {sequence_with_ends}") print() # Znajdź orientację CDS w pliku GenBank cds_orientation = None for feature in record.features: if feature.type == "CDS": cds_orientation = feature.location.strand break # Zakładamy, że interesuje nas pierwszy znaleziony CDS # Znajdź i wyświetl sekwencję promotora T7 z oznaczeniem na podstawie orientacji CDS for feature in record.features: if feature.type == "promoter" and "T7" in feature.qualifiers.get('note', [''])[0]: print("--") print("Promotor T7 znaleziony:") print(f"Type: {feature.type}") # Pobierz lokalizację start = feature.location.start end = feature.location.end # Pobierz sekwencję promotora promoter_seq = record.seq[start:end] # Oblicz sekwencję komplementarną complementary_seq = promoter_seq.complement() # Przygotuj format wyjścia na podstawie orientacji CDS if cds_orientation == 1: # CDS na nici dodatniej positive_strand_output = f" [{start}:{end}](+) 5'-{promoter_seq}-3'" negative_strand_output = f" [{start}:{end}](-) 3'-{complementary_seq}-5'" elif cds_orientation == -1: # CDS na nici ujemnej positive_strand_output = f" [{start}:{end}](+) 5'-{promoter_seq}-3'" negative_strand_output = f"* [{start}:{end}](-) 3'-{complementary_seq}-5'" else: # Orientacja CDS nieokreślona positive_strand_output = f" [{start}:{end}](unknown) 5'-{promoter_seq}-3'" negative_strand_output = f" [{start}:{end}](unknown) 3'-{complementary_seq}-5'" # Wyświetl wyniki print(f"(+): {positive_strand_output}") print(f"(-): {negative_strand_output}") print() #%% # Znajdź i wyświetl sekwencję terminatora T7 for feature in record.features: if feature.type == "terminator" and "T7" in feature.qualifiers.get('note', [''])[0]: print("Terminator T7 znaleziony:") print(f"Type: {feature.type}") print(f"Location: {feature.location}") # Pobierz sekwencję terminatora terminator_seq = record.seq[feature.location.start:feature.location.end] # Sprawdź orientację i odwróć/uzupełnij jeśli potrzebne if feature.strand == -1: terminator_seq = terminator_seq.reverse_complement() print(f"(+) Sequence: {terminator_seq}") print() #%% t=record[25:73].seq #%% a = str (t) a #%% class Model: complement = {'A': 'T', 'T': 'A', 'C': 'G', 'G': 'C'} def __init__ (self, seq=None): self.seq = seq def complement_dna(self): complementary_sequence = ''.join(complement[base] for base in self.seq) return complementary_sequence #%% r = Model(a) #%% r.seq #%% r.seq[::-1] #%% r.complement_dna() #%% r.seq #%% complement = {'A': 'T', 'T': 'A', 'C': 'G', 'G': 'C'} complementary_sequence = ''.join(complement[base] for base in r.seq) complementary_sequence #%% #%% # Sekwencja DNA dna_seq = record.seq # Znajdź miejsca cięcia dla NcoI i Bpu1102I ncoI_sites = NcoI.search(dna_seq) bpu1102I_sites = Bpu1102I.search(dna_seq) # Przyjmij pierwsze miejsca cięcia ncoI_site = ncoI_sites[0] bpu1102I_site = bpu1102I_sites[0] # Wyodrębnij fragmenty 20 zasad od miejsc cięcia ncoI_fragment = dna_seq[ncoI_site-16:ncoI_site+4] # 16 zasad przed i 4 po cięciu bpu1102I_fragment = dna_seq[bpu1102I_site:bpu1102I_site+20] # 20 zasad po cięciu # Lepkie końce ncoI_sticky_end = ncoI_fragment[-4:] ncoI_sticky_end_comp = ncoI_sticky_end.complement() bpu1102I_sticky_end = bpu1102I_fragment[:4] bpu1102I_sticky_end_comp = bpu1102I_sticky_end.complement() # Funkcja do wyświetlania dwuniciowego DNA z lepkimi końcami w formacie schodkowym z indeksami dla dłuższej nici def print_sticky_ends(seq1, sticky_end1, seq2, sticky_end2, start1, start2): # Convert sequences to strings seq1_str = str(seq1) sticky_end1_str = str(sticky_end1) comp_seq1_str = str(seq1.complement()) sticky_end1_comp_str = str(sticky_end1.complement()) seq2_str = str(seq2) sticky_end2_str = str(sticky_end2) comp_seq2_str = str(seq2.complement()) sticky_end2_comp_str = str(sticky_end2.complement()) # Add indexes to the ends seq1_with_index = f"{seq1_str[:-1]} -{start1+len(seq1_str)-1}" comp_seq1_with_index = f"{comp_seq1_str[:-1]}" seq2_with_index = f"{start2}- {seq2_str[1:]}" comp_seq2_with_index = f"{comp_seq2_str[1:]}" # Lustrzane odbicie sekwencji comp_seq2_with_index = comp_seq2_with_index[::-1] sticky_end2_comp_str = sticky_end2_comp_str[::-1] # Format the output as requested print(f"5'-{seq1_with_index} {sticky_end1_str}-3' (+)") print(f"3'-{comp_seq1_with_index}{sticky_end1_comp_str}-5'") print() print(f"5'-{sticky_end2_str} {seq2_with_index}-3' (+)") print(f" {sticky_end2_comp_str}{comp_seq2_with_index}-5'") # Wyświetl lepkie końce fragmentów print("Fragmenty po cięciu NcoI i Bpu1102I:") print_sticky_ends(ncoI_fragment[:-4], ncoI_sticky_end, bpu1102I_fragment[4:], bpu1102I_sticky_end, ncoI_site-16, bpu1102I_site) #%% # Projektowanie starterów forward_primer = gfp_seq[:20] # Pierwsze 20 nukleotydów sekwencji GFP reverse_primer = str(Seq(gfp_seq_with_his6[-20:]).reverse_complement()) # Ostatnie 20 nukleotydów + His6 print(f"Forward Primer: {forward_primer}") print(f"Reverse Primer: {reverse_primer}") #%% # Miejsca cięcia restryktazy enzymes = RestrictionBatch([EcoRI, BamHI, HindIII]) #enzymes = RestrictionBatch([BamHI, ]) restriction_sites = enzymes.search(record.seq) #%% # Przeprowadzenie analizy restrykcyjnej analysis = Analysis(AllEnzymes, record.seq) # Wynik analizy results = analysis.full() # Wyświetlanie wyników #for enzyme in results: # if len(results[enzyme]) > 0: # print(f"{enzyme}: {results[enzyme]}") #%% # Utwórz diagram plazmidu diagram = GenomeDiagram.Diagram("PUC19 Plasmid Map") track = diagram.new_track(1, name="Annotated Features", greytrack=True) feature_set = track.new_set() # Dodanie sekwencji kodującej GFP feature_set.add_feature( SeqFeature(FeatureLocation(0, len(gfp_seq)), strand=+1), name="GFP Coding Sequence", label=True, color=colors.lightblue ) # Dodanie His6 tagu feature_set.add_feature( SeqFeature(FeatureLocation(len(gfp_seq), len(gfp_seq_with_his6)), strand=+1), name="His6 Tag", label=True, color=colors.lightgreen ) # Dodanie forward primer feature_set.add_feature( SeqFeature(FeatureLocation(0, len(forward_primer)), strand=+1), name="Forward Primer", label=True, color=colors.orange ) # Dodanie reverse primer feature_set.add_feature( SeqFeature(FeatureLocation(len(gfp_seq_with_his6) - len(reverse_primer), len(gfp_seq_with_his6)), strand=-1), name="Reverse Primer", label=True, color=colors.red ) # Dodanie miejsc cięcia restryktazy for enzyme, sites in restriction_sites.items(): for site in sites: feature_set.add_feature( SeqFeature(FeatureLocation(site, site + 1), strand=0), name=enzyme, label=True, color=colors.purple ) # Rysowanie diagramu diagram.draw(format="circular", circular=True, pagesize='A4', start=0, end=len(record), circle_core=0.5) diagram.write("pdf/plasmid_map.pdf", "PDF") print("Zapisano diagram plazmidu do pliku 'plasmid_map.pdf'")