Кроим ДНК на Python — CRISPR gRNA finder, Часть I: Введение и базовый поиск

1. Часть I: Введение и базовый поиск (вы тута)

2. Часть II: Скоринг, off-target и работа с реальными данными (скоро)

Как я дошёл до жизни такой?

Несколько лет назад я наткнулся на статью про CRISPR-Cas9 ^[1] и домашние биолаборатории – люди буквально у себя дома экспериментировали с редактированием генов. Тогда это казалось чем-то из научной фантастики, но мысль засела. Потом появились новости о генной терапии ^[2], одобренные препараты (вот FDA ^[3], а у нас я только про BIOCAD слышал – СМА ^[4] и гемофилия B ^[5]), и тема перестала быть просто хайпом.

Я бэкенд-разработчик, книги по биологии в последний раз открывал в школе, но желание разобраться и посмотреть от этого не меньше. Проблема в том, что большинство материалов по биоинформатике либо требуют серьёзного биобэкграунда, либо сразу кидают в тебя терминами вроде “гомологичная рекомбинация” без объяснений.

Поэтому я решил зайти с другой стороны: взять конкретную задачу и разобраться в биологии ровно настолько, насколько нужно для её решения. А заодно написать инструмент, который можно будет потрогать руками.

Что будем делать?

Напишем CLI-утилиту для поиска потенциальных guide RNA ^[6] (gRNA) – это такие “наводчики” для CRISPR-Cas9, которые указывают, где именно резать ДНК.

Инструмент будет:

1. Читать последовательность ДНК из FASTA-файла ^[7]

2. Искать все возможные “сайты” (здесь и далее сайт – просто место, участок на молекуле ДНК/РНК) для CRISPR на обеих цепях.

3. Фильтровать кандидатов по GC-составу ^[8] и наличию повторов

4. Выводить результаты в удобном виде

Существуют готовые решения (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) ^[9].

Для самого популярного Cas9 (из бактерии Streptococcus pyogenes ^[10]) PAM – это NGG, где N – любой нуклеотид.

То есть структура выглядит так:

[20 нуклеотидов - сюда прилипает gRNA][NGG - PAM]

Нет PAM рядом – Cas9 там работать не будет.

Почему нельзя взять любые 20 букв?

Две причины:

1. PAM должен быть рядом – не везде есть последовательность NGG в нужном месте

2. 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

Логика ^[11] простая: ищем 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, без повторов.

Что дальше?

Статья получилась не очень большая, так как здесь я говорю только о базовых вещах, поэтому не обессудьте – обещаю, что в следующей части будет интереснее.

В следующей части мы:

1. Добавим скоринг кандидатов (не все gRNA с хорошим GC одинаково полезны)

2. Реализуем упрощённую проверку на off-target

3. Прогоним инструмент на реальном гене из NCBI

Код статьи ^[12]

---

P.S. Если вы биоинформатик и видите, что я где-то наврал – буду рад комментариям. Если вы разработчик и впервые слышите про CRISPR – надеюсь, было не слишком больно.

P.P.S Если есть желание предметнее обсудить тему статьи или обложить меня последними словами (по какой бы то ни было причине), то буду рад видеть вас тут: каналья ^[13].