#!/usr/bin/env python3 """ amplicon_digest.py Вырезает ампликон из полноразмерной последовательности гена по координатам праймеров (полученным из NCBI BLAST) и режет его указанной рестриктазой. Требуется: biopython pip install biopython --break-system-packages Как получить нужные параметры: 1. Сделайте BLAST праймеров против вашего gene.fna (или скачайте fasta с NCBI, где в заголовке указан диапазон геномных координат, например: >NC_000023.11:67544021-67730619 ...). 2. В выдаче BLAST найдите строки "Template" для forward и reverse праймера, например: Forward primer 1 CACAGGCTACCTGGTCCTGG 20 Template 67545458 .................... 67545477 Reverse primer 1 CTGCCTTACACAACTCCTTGGC 22 Template 67545873 ...................... 67545852 Здесь нужны 4 числа: fwd_start, fwd_end, rev_start, rev_end (координаты Template, в том порядке, как в выдаче). Пример запуска: python3 amplicon_digest.py \ --fasta gene.fna \ --fwd-start 67545458 --fwd-end 67545477 \ --rev-start 67545873 --rev-end 67545852 \ --enzyme StuI Если у вас несколько записей в fasta-файле (несколько сборок/контигов), укажите нужную через --record-index (0 - первая, по умолчанию 0) или --record-id (по подстроке в заголовке). Режим --custom (без BLAST-координат и fasta гена): Если у вас уже есть готовая последовательность ампликона (например, вы вручную подставили нужный аллель SNP), можно передать её напрямую: python3 amplicon_digest.py --custom ACGTACGT...GGCC --enzyme AspLEI Либо передать путь к fasta-файлу с уже готовым ампликоном: python3 amplicon_digest.py --custom amp2.fna --enzyme AspLEI В этом режиме --fasta и координаты праймеров не нужны. """ import argparse import sys from Bio import SeqIO from Bio.Restriction import Analysis, RestrictionBatch def parse_args(): p = argparse.ArgumentParser( description="Вырезать ампликон по координатам праймеров из BLAST " "и порезать его рестриктазой." ) p.add_argument("--fasta", default=None, help="Путь к fasta-файлу с полной последовательностью гена " "(не нужен, если используется --custom)") p.add_argument("--fwd-start", type=int, default=None, help="Координата начала forward-праймера (Template, из BLAST)") p.add_argument("--fwd-end", type=int, default=None, help="Координата конца forward-праймера (Template, из BLAST)") p.add_argument("--rev-start", type=int, default=None, help="Координата начала reverse-праймера (Template, из BLAST)") p.add_argument("--rev-end", type=int, default=None, help="Координата конца reverse-праймера (Template, из BLAST)") p.add_argument("--custom", default=None, help="Готовая последовательность ампликона, введённая прямо в консоли " "(строка ACGT...), ИЛИ путь к fasta-файлу с готовым ампликоном. " "Если указан этот флаг — --fasta и координаты праймеров не нужны, " "рестрикция считается сразу по этой последовательности. " "Удобно, например, чтобы вручную подставить SNP-аллель и " "посмотреть, режется сайт или нет.") p.add_argument("--enzyme", action="append", required=True, help="Название рестриктазы (например StuI). " "Можно указать несколько раз: --enzyme StuI --enzyme EcoRI") p.add_argument("--record-index", type=int, default=0, help="Индекс записи в fasta, если их несколько (по умолчанию 0)") p.add_argument("--record-id", default=None, help="Подстрока для поиска нужной записи по заголовку fasta " "(имеет приоритет над --record-index)") p.add_argument("--out-amplicon", default="amplicon.fasta", help="Куда сохранить fasta с ампликоном") return p.parse_args() def validate_args(args): if args.custom is None: missing = [name for name, val in [ ("--fasta", args.fasta), ("--fwd-start", args.fwd_start), ("--fwd-end", args.fwd_end), ("--rev-start", args.rev_start), ("--rev-end", args.rev_end), ] if val is None] if missing: sys.exit("Ошибка: без --custom нужно указать все из: " "--fasta, --fwd-start, --fwd-end, --rev-start, --rev-end. " f"Не хватает: {', '.join(missing)}") def load_custom_sequence(custom_arg): """ --custom может быть либо путём к fasta-файлу с готовым ампликоном, либо просто "сырой" строкой ACGT..., введённой прямо в консоли. """ import os from Bio.Seq import Seq if os.path.isfile(custom_arg): record = SeqIO.read(custom_arg, "fasta") print(f"Загружен готовый ампликон из файла: {custom_arg} " f"(id={record.id}, длина {len(record.seq)} п.н.)") return record.seq raw = "".join(custom_arg.split()).upper() # убираем пробелы/переносы строк invalid_chars = set(raw) - set("ACGTNRYSWKMBDHV") if invalid_chars: sys.exit(f"Ошибка: в --custom встречены недопустимые символы: {invalid_chars}. " "Ожидается последовательность нуклеотидов (ACGT...).") print(f"Используется кастомная последовательность, введённая в консоли " f"(длина {len(raw)} п.н.)") return Seq(raw) def parse_header_coords(header: str): """ Пытается вытащить геномные координаты и направление из заголовка вида: 'NC_000023.11:67544021-67730619 ...' (плюс-цепь, координаты растут) 'NC_000010.11:c121598458-121478330 ...' (минус-цепь, префикс 'c' = complement, координаты в заголовке убывают) Возвращает (start_coord, is_complement) или (None, None), если не удалось распарсить. start_coord — геномная координата, соответствующая ПЕРВОЙ позиции (файл-позиция 1) последовательности в fasta-файле. is_complement — True, если последовательность идёт по убыванию геномных координат (т.е. в файле уже лежит обратный комплемент участка генома). """ try: coords_part = header.split()[0] # NC_000010.11:c121598458-121478330 coords_part = coords_part.split(":")[-1] # c121598458-121478330 is_complement = coords_part.startswith("c") if is_complement: coords_part = coords_part[1:] # 121598458-121478330 start_str = coords_part.split("-")[0] return int(start_str), is_complement except (IndexError, ValueError): return None, None def load_record(fasta_path, record_index, record_id): records = list(SeqIO.parse(fasta_path, "fasta")) if not records: sys.exit(f"Ошибка: в файле {fasta_path} не найдено ни одной fasta-записи.") if record_id: for r in records: if record_id in r.description: return r sys.exit(f"Ошибка: запись с '{record_id}' в заголовке не найдена.") if record_index >= len(records): sys.exit(f"Ошибка: в файле только {len(records)} записей, " f"а запрошен индекс {record_index}.") return records[record_index] def main(): args = parse_args() validate_args(args) if args.custom is not None: amplicon = load_custom_sequence(args.custom) print(f"\nПоследовательность:\n{amplicon}") with open(args.out_amplicon, "w") as f: f.write(">custom_amplicon\n") f.write(str(amplicon) + "\n") print(f"\nСохранено в: {args.out_amplicon}") else: record = load_record(args.fasta, args.record_index, args.record_id) seq = record.seq print(f"Загружена последовательность: {record.id} (длина {len(seq)} п.н.)") gene_start_genomic, is_complement = parse_header_coords(record.description) if gene_start_genomic is None: print("Не удалось определить геномную стартовую координату из заголовка.") print("Предполагаю, что переданные координаты уже соответствуют " "позициям внутри файла (нумерация с 1).") gene_start_genomic = 1 is_complement = False else: strand = "минус (complement, 'c' в заголовке)" if is_complement else "плюс" print(f"Геномная стартовая координата последовательности: {gene_start_genomic} " f"(цепь: {strand})") def g2p(g): """Геномная координата -> 1-based позиция внутри файла.""" if is_complement: # координаты в заголовке убывают: файл-позиция 1 = gene_start_genomic, # дальше по файлу геномная координата уменьшается return gene_start_genomic - g + 1 return g - gene_start_genomic + 1 fwd_start = g2p(args.fwd_start) fwd_end = g2p(args.fwd_end) rev_start = g2p(args.rev_start) rev_end = g2p(args.rev_end) # Праймеры могут быть на любой цепи, поэтому берём min/max для границ ампликона amplicon_left = min(fwd_start, fwd_end, rev_start, rev_end) amplicon_right = max(fwd_start, fwd_end, rev_start, rev_end) if amplicon_left < 1 or amplicon_right > len(seq): sys.exit("Ошибка: вычисленные координаты ампликона выходят за пределы " "последовательности. Проверьте правильность введённых координат " "и совпадение записи fasta с той, что использовалась в BLAST.") amplicon = seq[amplicon_left - 1:amplicon_right] print(f"\nКоордината ампликона внутри файла: {amplicon_left}-{amplicon_right}") print(f"Длина ампликона: {len(amplicon)} п.н.") print(f"\nПоследовательность ампликона:\n{amplicon}") with open(args.out_amplicon, "w") as f: f.write(f">amplicon_{record.id}_{args.fwd_start}-{args.rev_start}\n") f.write(str(amplicon) + "\n") print(f"\nАмпликон сохранён в: {args.out_amplicon}") # --- Рестрикция --- from Bio.Restriction import AllEnzymes all_names = {str(e): str(e) for e in AllEnzymes} all_names_lower = {name.lower(): name for name in all_names} resolved_enzymes = [] for requested in args.enzyme: if requested in all_names: resolved_enzymes.append(requested) elif requested.lower() in all_names_lower: correct_name = all_names_lower[requested.lower()] print(f"Внимание: фермент '{requested}' не найден, но нашлось " f"совпадение без учёта регистра: '{correct_name}'. Использую его.") resolved_enzymes.append(correct_name) else: # ищем похожие названия, чтобы подсказать candidates = [n for n in all_names if requested.lower() in n.lower() or n.lower() in requested.lower()] msg = f"Ошибка: фермент '{requested}' не найден в базе REBASE (Biopython)." if candidates: msg += f" Возможно, вы имели в виду: {', '.join(candidates[:10])}" sys.exit(msg) rb = RestrictionBatch(resolved_enzymes) ana = Analysis(rb, amplicon) print("\n=== Сайты рестрикции ===") ana.print_that() for enzyme in rb: fragments = enzyme.catalyze(amplicon) print(f"\n=== Фрагменты после {enzyme} ===") if len(fragments) == 1: print("Сайтов не найдено, ампликон не режется.") continue for i, frag in enumerate(fragments, 1): print(f"Фрагмент {i}: {len(frag)} п.н.") print(str(frag)) if __name__ == "__main__": main()