pfas/py/protocol.py

#! /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'")