506 lines
17 KiB
Python
506 lines
17 KiB
Python
#! /usr/bin/env python3
|
|
# vim:fenc=utf-8
|
|
#
|
|
# Copyright © 2024 user <mpabi@mpabi.pl>
|
|
#
|
|
# 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'")
|