bioinf/py/3.py

69 lines
2.6 KiB
Python
Raw Normal View History

2024-10-23 17:08:00 +00:00
from Bio import Entrez, SeqIO
from Bio.Seq import Seq
# Ustaw swój adres e-mail (wymagany przez NCBI)
Entrez.email = "your_email@example.com"
# Funkcja do pobierania szczegółów genu
def get_gene_info(gene_id):
try:
# Pobierz adnotacje genowe z NCBI
handle = Entrez.efetch(db="gene", id=gene_id, retmode="xml")
records = Entrez.read(handle)
# Wyciągnij podstawowe informacje
gene_name = records[0]["Entrezgene_gene"]["Gene-ref"]["Gene-ref_locus"]
summary = records[0]["Entrezgene_summary"]
chromosome = records[0]["Entrezgene_source"]["BioSource"]["BioSource_subtype"][0]["SubSource_name"]
exon_count = len(records[0]["Entrezgene_locus"][0]["Gene-commentary_seqs"])
# Wyświetl podstawowe informacje o genie
print(f"Nazwa genu: {gene_name}")
print(f"Opis: {summary}")
print(f"Chromosom: {chromosome}")
print(f"Liczba eksonów: {exon_count}")
except Exception as e:
print(f"Wystąpił błąd podczas pobierania informacji o genie: {e}")
# Funkcja do załadowania sekwencji genów z pliku FASTA
def load_sequences(file_path):
return list(SeqIO.parse(file_path, "fasta"))
# Znajdź potencjalne sekwencje gRNA (20 nukleotydów + PAM "NGG") w określonym regionie
def find_crispr_sites_in_region(sequence, start, end):
pam = "NGG"
sites = []
# Ogranicz przeszukiwanie do wybranego regionu
target_region = sequence[start:end]
for i in range(len(target_region) - 23):
target = target_region[i:i+20]
pam_site = target_region[i+20:i+23]
if pam_site.endswith("GG"):
# Zapisz pełną sekwencję wraz z lokalizacją w genomie
sites.append((start + i, target + pam_site))
return sites
# Przykład użycia
gene_id = "835067" # Przykładowe ID genu (należy zastąpić rzeczywistym ID odpornego na suszę genu)
get_gene_info(gene_id)
# Zakładamy, że sekwencja DNA jest już wczytana do pliku FASTA
genes = load_sequences("arabidopsis_genome.fasta")
# Przykładowe pozycje startowe i końcowe genu odpowiedzialnego za odporność na suszę
start_position = 1000 # przykładowa pozycja początkowa
end_position = 2000 # przykładowa pozycja końcowa
for gene in genes:
print(f"\nGene ID: {gene.id}")
crispr_sites = find_crispr_sites_in_region(str(gene.seq), start_position, end_position)
# Wyświetl tylko 5 pierwszych wyników, wraz z ich pozycją w sekwencji DNA
for pos, site in crispr_sites[:5]:
print(f"Position: {pos}, CRISPR site: {site}")