Кроим ДНК на Python — CRISPR gRNA finder, Часть I: Введение и базовый поиск. cas9.. cas9. crispr.. cas9. crispr. анализ днк.. cas9. crispr. анализ днк. биоинженерия.. cas9. crispr. анализ днк. биоинженерия. биоинформатика.. cas9. crispr. анализ днк. биоинженерия. биоинформатика. днк.. cas9. crispr. анализ днк. биоинженерия. биоинформатика. днк. чайник-чайнику.

CRISPR-Cas9

CRISPR-Cas9

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

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

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

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

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

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

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

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

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

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

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

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

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

Для самого популярного Cas9 (из бактерии Streptococcus pyogenes) 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

Логика простая: ищем 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

Код статьи

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

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

Автор: egnodus

Источник

Версия для печати Версия для печати
Rambler's Top100