Часть I: Введение и базовый поиск (вы тута)
Как я дошёл до жизни такой?
Несколько лет назад я наткнулся на статью про CRISPR-Cas9 и домашние биолаборатории - люди буквально у себя дома экспериментировали с редактированием генов. Тогда это казалось чем-то из научной фантастики, но мысль засела. Потом появились новости о генной терапии, одобренные препараты (вот FDA, а у нас я только про BIOCAD слышал - СМА и гемофилия B), и тема перестала быть просто хайпом.
Я бэкенд-разработчик, книги по биологии в последний раз открывал в школе, но желание разобраться и посмотреть от этого не меньше. Проблема в том, что большинство материалов по биоинформатике либо требуют серьёзного биобэкграунда, либо сразу кидают в тебя терминами вроде "гомологичная рекомбинация" без объяснений.
Поэтому я решил зайти с другой стороны: взять конкретную задачу и разобраться в биологии ровно настолько, насколько нужно для её решения. А заодно написать инструмент, который можно будет потрогать руками.
Что будем делать?
Напишем CLI-утилиту для поиска потенциальных guide RNA (gRNA) - это такие "наводчики" для CRISPR-Cas9, которые указывают, где именно резать ДНК.
Инструмент будет:
Читать последовательность ДНК из FASTA-файла
Искать все возможные "сайты" (здесь и далее сайт - просто место, участок на молекуле ДНК/РНК) для CRISPR на обеих цепях.
Фильтровать кандидатов по GC-составу и наличию повторов
Выводить результаты в удобном виде
Существуют готов��е решения (CRISPOR, Benchling, CRISPRdirect), но наша задача просто понять, как они работают изнутри, а не конкурировать с ними.
P.S. Если вы, как и я, последний раз слышали про ДНК на уроке биологии - не переживайте, разберёмся вместе. Обещаю минимум терминов и максимум кода.
Минимум биологии (правда минимум)
Чтобы писать код, нужно понять несколько базовых вещей. Постараюсь объяснить без занудства.
ДНК - это строчка из 4 букв
ДНК - две цепочки, скрученные в спираль. Каждая цепочка - последовательность из четырёх нуклеотидов: A (аденин), T (тимин), G (гуанин), C (цитозин).
Цепочки комплементарны друг другу, то есть значение в одной цепочке связано со значением в другой на той же позиции, как то: A всегда напротив T, G всегда напротив C. Если одна цепочка AATGC, вторая - TTACG.
Когда записывают последовательность гена, обычно пишут только одну цепочку. Вторую можно вычислить через словарь подстановки (как в шифре Цезаря).
CRISPR-Cas9 - молекулярные ножницы
CRISPR-Cas9 состоит из двух частей:
Cas9 - белок, который физически разрезает обе цепи ДНК
guide RNA (gRNA) - короткая последовательность (~20 нуклеотидов), которая говорит Cas9, куда идти
gRNA комплементарна целевому участку ДНК. Она "прилипает" к нужному месту, Cas9 режет рядом.
PAM - обязательное условие
Cas9 не может резать где угодно. Ему нужно, чтобы сразу после целевого участка была специальная последовательность - PAM (Protospacer Adjacent Motif).
Для самого популярного Cas9 (из бактерии Streptococcus pyogenes) PAM - это NGG, где N - любой нуклеотид.
То есть структура выглядит так:
[20 нуклеотидов - сюда прилипает gRNA][NGG - PAM]
Нет PAM рядом - Cas9 там работать не будет.
Почему нельзя взять любые 20 букв?
Две причины:
PAM должен быть рядом - не везде есть последовательность NGG в нужном месте
Off-target эффекты - геном человека содержит 3 миллиарда пар нуклеотидов. Последовательность из 20 букв может встретиться не один раз. Если gRNA "прилипнет" не туда - Cas9 разрежет не тот ген
Поэтому при дизайне gRNA важно проверять уникальность последовательности.
Критерии хорошей gRNA
Кроме наличия PAM, есть ещё несколько критериев:
GC-состав - соотношение G и C в последовательности. Оптимально 40-60%. Слишком низкий - gRNA плохо связывается. Слишком высокий - может связываться неспецифично.
Повторы - если в gRNA есть AAAA или TTTT - это плохо. РНК-полимераза может "соскользнуть", gRNA будет синтезироваться с ошибками.
Off-target - нужно проверить, не встречается ли последовательность где-то ещё в геноме. Это самый важный критерий, но и самый сложный - требует поиска по всему геному. В этой части мы его пропустим, вернёмся во второй.
Пишем код
Хватит теории, погнали кодить.
Поиск PAM-сайтов
Начнём с функции, которая находит все позиции PAM (NGG) в последовательности:
def find_ngg_pam(sequence: str) -> list[int]: """Находит все позиции PAM (NGG) в последовательности.""" candidates_positions = [] last_idx = -1 while True: last_idx = sequence.find('GG', last_idx + 1) if last_idx > 0: candidates_positions.append(last_idx - 1) elif last_idx == -1: break return candidates_positions
Логика простая: ищем GG и берём позицию на один символ раньше (там N, то есть любой из A, C, T, G). Проверка last_idx > 0 отсекает случай, когда GG в самом начале строки, так как перед ним нет нуклеотида.
Извлечение кандидатов на роль gRNA
Теперь для каждого PAM берём 20 нуклеотидов перед ним:
def extract_grna_candidates(sequence: str, pam_positions: list[int]) -> list[dict]: """Извлекает кандидатов на gRNA (20 нуклеотидов перед PAM).""" grna_candidates = [] for pam_pos in pam_positions: if pam_pos >= 20: grna_seq = sequence[pam_pos - 20:pam_pos] grna_candidates.append({ 'sequence': grna_seq, 'pam_position': pam_pos }) return grna_candidates
Если перед PAM меньше 20 нуклеотидов - пропускаем, там не получится полноценная gRNA.
GC-состав
def calculate_gc_content(sequence: str) -> float: """Считает процент G и C в последовательности.""" gc_count = sequence.count('G') + sequence.count('C') return (gc_count / len(sequence)) * 100 if len(sequence) > 0 else 0.0
Поиск повторов
def has_poly_repeats(sequence: str, min_length: int = 4) -> bool: """Проверяет наличие повторов одного нуклеотида.""" for base in 'ACGT': if base * min_length in sequence: return True return False
Reverse complement
ДНК двухцепочечная, CRISPR может работать на обеих цепях. Чтобы искать на второй цепи, нужно сделать reverse complement - заменить буквы на комплементарные и перевернуть строку:
def reverse_complement(sequence: str) -> str: """Возвращает обратную комплементарную последовательность.""" complement = str.maketrans('ACGT', 'TGCA') return sequence.translate(complement)[::-1]
Собираем всё вместе
def _analyze_sequence(sequence: str, is_forward: bool) -> list[dict]: """Анализирует одну цепь.""" pam_positions = find_ngg_pam(sequence) grna_candidates = extract_grna_candidates(sequence, pam_positions) analyzed_candidates = [] for candidate in grna_candidates: gc_content = calculate_gc_content(candidate['sequence']) poly_repeats = has_poly_repeats(candidate['sequence']) analyzed_candidates.append({ 'sequence': candidate['sequence'], 'pam_position': candidate['pam_position'], 'gc_content': gc_content, 'has_poly_repeats': poly_repeats, 'is_forward': is_forward }) return analyzed_candidates def analyze_sequence(sequence: str) -> list[dict]: """Анализирует обе цепи ДНК.""" candidates = [] candidates.extend(_analyze_sequence(sequence, is_forward=True)) rev_comp_sequence = reverse_complement(sequence) candidates.extend(_analyze_sequence(rev_comp_sequence, is_forward=False)) return candidates
Фильтрация
def filter_candidates( candidates: list[dict], gc_min: float = 40.0, gc_max: float = 60.0, allow_repeats: bool = False ) -> list[dict]: """Фильтрует кандидатов по критериям качества.""" filtered = [] for candidate in candidates: if gc_min <= candidate['gc_content'] <= gc_max: if allow_repeats or not candidate['has_poly_repeats']: filtered.append(candidate) return filtered
CLI-интерфейс
Обернём всё в CLI с помощью typer. Заодно добавим парсер FASTA-файлов:
import typer from pathlib import Path app = typer.Typer() def parse_fasta(file_path: Path) -> dict[str, str]: """Парсит FASTA файл, возвращает {название: последовательность}.""" if not file_path.exists(): typer.echo(f"Файл не найден: {file_path}") raise typer.Exit(code=1) raw_data = file_path.read_text() results = {} sequence = "" name = "" for line in raw_data.splitlines(): if line.startswith(">"): if sequence and name: results[name] = sequence name = line[1:].strip() sequence = "" else: sequence += line.strip().upper() if sequence and name: results[name] = sequence return results @app.command() def find( fasta_path: Path = typer.Argument(..., help="Путь к FASTA файлу"), gc_min: float = typer.Option(40.0, help="Минимальный GC-состав"), gc_max: float = typer.Option(60.0, help="Максимальный GC-состав"), allow_repeats: bool = typer.Option(False, help="Разрешить повторы"), ): """Поиск CRISPR gRNA кандидатов в последовательности.""" sequences = parse_fasta(fasta_path) for name, sequence in sequences.items(): typer.echo(f"\n=== {name} ===") candidates = analyze_sequence(sequence) filtered = filter_candidates(candidates, gc_min, gc_max, allow_repeats) if not filtered: typer.echo("Кандидатов не найдено") continue for c in filtered: strand = "+" if c["is_forward"] else "-" typer.echo( f" {c['sequence']} | pos: {c['pam_position']} | " f"GC: {c['gc_content']:.1f}% | strand: {strand}" ) if __name__ == "__main__": app()
Проверяем
Создадим тестовый FASTA-файл:
>test_sequence GCCTACGGATCCTGAAGGCTTAGCCTGAGGAATCG
Запускаем:
python crispr_finder.py test.fasta
Получаем:
=== test_sequence === GATCCTGAAGGCTTAGCCTG | pos: 27 | GC: 55.0% | strand: + ATTCCTCAGGCTAAGCCTTC | pos: 22 | GC: 50.0% | strand: - GCTAAGCCTTCAGGATCCGT | pos: 31 | GC: 55.0% | strand: -
Один кандидат на прямой цепи, два на обратной. Все в диапазоне 40-60% GC, без повторов.
Что дальше?
Статья получилась не очень большая, так как здесь я говорю только о базовых вещах, поэтому не обессудьте - обещаю, что в следующей части будет интереснее.
В следующей части мы:
Добавим скоринг кандидатов (не все gRNA с хорошим GC одинаково полезны)
Реализуем упрощённую проверку на off-target
Прогоним инструмент на реальном гене из NCBI
P.S. Если вы биоинформатик и видите, что я где-то наврал - буду рад комментариям. Если вы разработчик и впервые слышите про CRISPR - надеюсь, было не слишком больно.
P.P.S Если есть желание предметнее обсудить тему статьи или обложить меня последними словами (по какой бы то ни было причине), то буду рад видеть вас тут: каналья.
