Часть II: Скоринг, off-target и работа с реальными данными (вы тута)
В первой части мы разобрались с базовой биологией CRISPR-Cas9 и написали CLI-утилиту для поиска gRNA. Фильтровали по GC-составу и наличию повторов. В итоге получился рабочий, но очень простой инструмент.
Сегодня добавим функций: научим программу сортировать кандидатов по качеству и проверять, не порежет ли наша gRNA что-нибудь лишнее в геноме. А чтобы было интереснее прогоним всё на настоящем гене с криминальной историей.
Почему CCR5?
Помните скандал 2018 года? Китайский учёный Хэ Цзянькуй отредактировал геном человеческих эмбрионов и довёл беременность до родов. На свет появились близнецы Лулу и Нана - первые генетически модифицированные люди. Учёный получил три года тюрьмы, а научное сообщество - этическую мигрень на десятилетие вперёд.
Хэ целился в ген CCR5. Этот ген кодирует белок-рецептор на поверхности иммунных клеток - это что-то вроде дверной ручки, за которую ВИЧ хватается, чтобы проникнуть внутрь. Примерно у 1% европейцев есть мутация CCR5-Δ32: у них эта "ручка" сломана и, как следствие, такие люди устойчивы к ВИЧ от рождения.
Идея Хэ была в том, чтобы искусственно "сломать" CCR5 у эмбрионов и сделать детей невосприимчивыми к вирусу. Этика эксперимента - это отдельный разговор, но ген CCR5 с тех пор стал, пожалуй, самой известной мишенью для CRISPR.
Его и возьмём для тестов.
Ликбез: 5' и 3' - что за штрихи?
Прежде чем лезть в код, разберёмся с одной штукой, которая постоянно мелькает в биоинформатике.
ДНК - это цепочка нуклеотидов (A, C, T, G), связанных через сахарофосфатный "скелет". Каждый сахар в этом скелете - молекула дезоксирибозы с пронумерованными атомами углерода: 1', 2', 3', 4', 5' (читается "один штрих", "пять штрих"). Штрих - чтобы не путать с номерами атомов в азотистых основаниях.
Цепь ДНК всегда имеет направление, как улица с односторонним движением:
5'-конец - начало цепи, там торчит свободная фосфатная группа
3'-конец - конец цепи, там свободная OH-группа
Представьте цепочку вагонов: 5'-конец - это локомотив, 3'-конец - последний вагон. Все ферменты "читают" ДНК в одном направлении (от 5' к 3'), как поезд едет от головы к хвосту.
Две цепи ДНК идут "навстречу" друг другу:
5'-ATGCCTGA-3' → (прямая цепь)
||||||||
3'-TACGGACT-5' ← (комплементарная)
Когда мы говорим "3'-конец gRNA", имеем в виду ту часть, которая ближе к PAM. Это важно для скоринга - скоро увидите почему.
Скоринг: не все gRNA одинаково полезны
В первой части мы отбирали gRNA по двум критериям: GC от 40 до 60% и отсутствие повторов. Но внутри этого диапазона кандидаты сильно различаются по эффективности.
Представьте, что вы выбираете сотрудника. Формально все кандидаты подходят (есть диплом, опыт работы), но кто-то явно лучше, поэтому нужна система оценки.
Что влияет на качество gRNA
GC-состав около 50% - золотая середина. Слишком много GC - цепи слипаются слишком крепко, gRNA не хочет отпускать ДНК после разрезания. Слишком мало - связывание слабое, Cas9 может "соскочить" раньше времени.
Нуклеотиды на концах - 3'-конец gRNA (ближе к PAM) критичен для связывания. Если там много G и C, комплекс получается слишком стабильным - как клей, который схватился намертво. Cas9 режет, но не может отцепиться и пойти дальше.
Первый нуклеотид - gRNA обычно экспрессируют с U6-промотора, а он капризный. Если первая буква - T (в РНК это U), промотор работает плохо. Если G - отлично.
Формула скоринга
Начинаем со 100 очков и вычитаем штрафы:
def calculate_score(sequence: str, gc_content: float, has_repeats: bool) -> float:
"""
Оценивает качество gRNA.
100 = идеал, штрафы уменьшают скор.
Бонус за G в начале может дать >100.
"""
score = 100.0
# Штраф за отклонение от оптимального GC (50%)
# Каждый процент отклонения = -1.5 очка
score -= abs(gc_content - 50) * 1.5
# повторы - это плохо
if has_repeats:
score -= 20
# Проверяем 3'-конец (последние 6 нуклеотидов)
last_6 = sequence[-6:]
gc_at_end = last_6.count('G') + last_6.count('C')
if gc_at_end > 3:
score -= (gc_at_end - 3) * 5
# Первый нуклеотид
if sequence[0] == 'T':
score -= 10 # U6 не любит тимин
elif sequence[0] == 'G':
score += 5 # U6 любит гуанин
return max(0, score)
Это упрощённая модель. Настоящие алгоритмы вроде Rule Set 2 используют машинное обучение на тысячах экспериментов, учитывают позицию каждого нуклеотида и ещё кучу факторов. Но для понимания принципа - достаточно.
Off-target: ищем случайных жертв
Самая большая проблема CRISPR специфичность. Геном человека - это три миллиарда пар нуклеотидов. Наша gRNA составляет всего 20 букв. Какова вероятность, что такая последовательность встретится где-то ещё?
Спойлер: высокая. И даже если точного совпадения нет, Cas9 может разрезать участки с 1-2 отличиями. Представьте, что вы ищете в толпе человека по описанию "высокий брюнет в синей куртке". Найдёте нужного, но по дороге обратите внимание на пару похожих.
Для поиска таких "похожих" используют BLAST.
Что такое BLAST
BLAST (Basic Local Alignment Search Tool) - это Google для биологов. Вы даёте ему последовательность, он ищет похожие в гигантских базах данных: геномах, транскриптомах, протеомах.
Алгоритм хитрый: вместо того чтобы сравнивать каждую букву (это заняло бы вечность), BLAST сначала ищет короткие "слова" такие совпадения из 7-11 букв. Нашёл слово и расширяет в обе стороны, проверяя, насколько длинное совпадение получится.
В ответ отдаёт список отобранных и метрики:
E-value - вероятность случайного совпадения. E-value = 0.001 значит "такое совпадение случайно встретится раз на 1000 запросов". Чем меньше - тем надёжнее.
Identity - процент идентичных букв. 100% = точное совпадение.
Coverage - какая часть запроса покрыта совпадением.
Biopython: швейцарский нож биоинформатика
Biopython - библиотека, которая умеет почти всё: читать форматы последовательностей, работать с базами NCBI, запускать BLAST, строить филогенетические деревья и многое другое.
Нам нужны два модуля:
Bio.Blast.NCBIWWW- отправляет запросы на сервер NCBIBio.Blast.NCBIXML- парсит результаты
pip install biopython
Проверка off-target через BLAST
from Bio.Blast import NCBIWWW, NCBIXML
def check_off_targets(
sequence: str,
organism: str = "Homo sapiens",
max_hits: int = 50,
) -> list[dict]:
"""
Ищет потенциальные off-target через NCBI BLAST.
Внимание: запрос занимает 30-60 секунд!
NCBI просит не частить - максимум 1 запрос в 10 секунд.
"""
# Добавляем PAM к запросу
query = sequence + "NGG"
result_handle = NCBIWWW.qblast(
program="blastn", # нуклеотидный поиск
database="nt", # nucleotide collection
sequence=query,
entrez_query=f'"{organism}"[organism]',
hitlist_size=max_hits,
word_size=7, # для коротких последовательностей
expect=1000, # высокий порог, чтобы поймать слабые совпадения
megablast=False
)
off_targets = []
for record in NCBIXML.parse(result_handle):
for alignment in record.alignments:
for hsp in alignment.hsps:
identity = (hsp.identities / hsp.align_length) * 100
off_targets.append({
'title': alignment.title[:80],
'identity': round(identity, 1),
'mismatches': hsp.align_length - hsp.identities,
'e_value': hsp.expect,
})
result_handle.close()
return off_targets
Пара замечаний:
word_size=7 - обычно BLAST ищет "слова" из 11 букв, но наш запрос всего 23 буквы. С большим word_size он может ничего не найти.
expect=1000 - нарочно высокий порог. Для off-target нас интересуют даже слабые совпадения: участок с 85% идентичности всё ещё может быть разрезан Cas9.
Оценка риска
Что считать опасным?
≥95% идентичности - Cas9 скорее всего разрежет
85-95% - может разрезать, зависит от позиции мисматчей
<85% - обычно безопасно
def assess_off_target_risk(
off_targets: list[dict],
target_gene: str
) -> dict:
"""Сортирует off-target по уровню риска."""
# Паттерны для фильтрации целевого гена
short_name = target_gene.split()[0].split('_')[0].lower()
target_patterns = [
short_name,
"chemokine receptor 5",
"cc chemokine receptor",
"c-c motif chemokine receptor",
"ccr-5",
"nonfunctional cc chemokine",
"cc chemokine re",
]
high_risk = []
medium_risk = []
low_risk = []
for ot in off_targets:
title_lower = ot['title'].lower()
# Пропускаем сам целевой ген - это не off-target
if any(p in title_lower for p in target_patterns):
continue
if ot['identity'] >= 95:
high_risk.append(ot)
elif ot['identity'] >= 85:
medium_risk.append(ot)
else:
low_risk.append(ot)
return {
'high_risk': high_risk,
'medium_risk': medium_risk,
'low_risk': low_risk,
'is_safe': len(high_risk) == 0,
}
Да-да, знаю, что target_patterns убивает переиспользование для других генов, не CCR5, но для примера никакого простого способа собрать все синонимы, по имени из FASTA файла, я не нашёл :(
Собираем всё вместе
Код второй части живёт в отдельном файле и импортирует функции из первой:
# crispr_finder_v2.py
from pathlib import Path
from Bio.Blast import NCBIWWW, NCBIXML
import typer
from crispr_finder import (
parse_fasta,
analyze_sequence,
filter_candidates,
)
app = typer.Typer()
def calculate_score(sequence: str, gc_content: float, has_repeats: bool) -> float:
"""Оценивает качество gRNA по шкале 0-100+."""
score = 100.0
score -= abs(gc_content - 50) * 1.5
if has_repeats:
score -= 20
last_6 = sequence[-6:]
gc_at_end = last_6.count('G') + last_6.count('C')
if gc_at_end > 3:
score -= (gc_at_end - 3) * 5
if sequence[0] == 'T':
score -= 10
elif sequence[0] == 'G':
score += 5
return max(0, score)
def score_candidates(candidates: list[dict]) -> list[dict]:
"""Добавляет скор к каждому кандидату."""
for c in candidates:
c['score'] = calculate_score(
c['sequence'],
c['gc_content'],
c['has_poly_repeats']
)
return candidates
def check_off_targets(
sequence: str,
organism: str = "Homo sapiens",
max_hits: int = 50,
) -> list[dict]:
"""Ищет off-target через NCBI BLAST."""
query = sequence + "NGG"
result_handle = NCBIWWW.qblast(
program="blastn",
database="nt",
sequence=query,
entrez_query=f'"{organism}"[organism]',
hitlist_size=max_hits,
word_size=7,
expect=1000,
megablast=False
)
off_targets = []
for record in NCBIXML.parse(result_handle):
for alignment in record.alignments:
for hsp in alignment.hsps:
identity = (hsp.identities / hsp.align_length) * 100
off_targets.append({
'title': alignment.title[:80],
'identity': round(identity, 1),
'mismatches': hsp.align_length - hsp.identities,
'e_value': hsp.expect,
})
result_handle.close()
return off_targets
def assess_off_target_risk(off_targets: list[dict], target_gene: str) -> dict:
"""Оценивает риск off-target эффектов."""
short_name = target_gene.split()[0].split('_')[0].lower()
target_patterns = [
short_name,
"chemokine receptor 5",
"c-c motif chemokine receptor 5",
]
high_risk = []
medium_risk = []
low_risk = []
for ot in off_targets:
title_lower = ot['title'].lower()
if any(p in title_lower for p in target_patterns):
continue
if ot['identity'] >= 95:
high_risk.append(ot)
elif ot['identity'] >= 85:
medium_risk.append(ot)
else:
low_risk.append(ot)
return {
'high_risk': high_risk,
'medium_risk': medium_risk,
'low_risk': low_risk,
'is_safe': len(high_risk) == 0,
}
@app.command()
def analyze(
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-состав"),
top_n: int = typer.Option(10, help="Показать топ N кандидатов"),
check_offtarget: bool = typer.Option(False, help="Проверить off-target (медленно!)"),
):
"""Анализ gRNA с скорингом и опциональной проверкой off-target."""
sequences = parse_fasta(fasta_path)
for name, sequence in sequences.items():
typer.echo(f"\n{'=' * 60}")
typer.echo(f"Ген: {name}")
typer.echo(f"Длина: {len(sequence)} нуклеотидов")
typer.echo(f"{'=' * 60}")
# Ищем и фильтруем кандидатов (функции из первой части)
candidates = analyze_sequence(sequence)
filtered = filter_candidates(candidates, gc_min, gc_max, allow_repeats=False)
# Добавляем скоринг
scored = score_candidates(filtered)
scored.sort(key=lambda x: x['score'], reverse=True)
if not scored:
typer.echo("Кандидатов не найдено")
continue
typer.echo(f"\nНайдено кандидатов: {len(scored)}")
typer.echo(f"Топ-{min(top_n, len(scored))} по скору:\n")
for i, c in enumerate(scored[:top_n], 1):
strand = "+" if c['is_strand'] else "-"
typer.echo(
f"{i:2}. {c['sequence']} | "
f"score: {c['score']:.1f} | "
f"GC: {c['gc_content']:.1f}% | "
f"strand: {strand}"
)
if check_offtarget:
typer.echo(" Проверяю off-target...")
try:
off_targets = check_off_targets(c['sequence'])
risk = assess_off_target_risk(off_targets, name)
if risk['is_safe']:
typer.echo(" ✓ Высокий риск off-target не обнаружен")
else:
typer.echo(
f" ⚠ Риск: {len(risk['high_risk'])} высокий, "
f"{len(risk['medium_risk'])} средний"
)
for ot in risk['high_risk'][:3]:
typer.echo(f" - {ot['title']}: {ot['identity']}%")
except Exception as e:
typer.echo(f" ✗ Ошибка: {e}")
if __name__ == "__main__":
app()
Тестируем на CCR5
Скачиваем последовательность CCR5 (или берём из репозитория) и запускаем:
python crispr_finder_v2.py analyze ccr5.fasta --top-n 5
============================================================
Ген: CCR5_CDS Homo sapiens C-C motif chemokine receptor 5 (CCR5), mRNA
Длина: 1059 нуклеотидов
============================================================
Найдено кандидатов: 62
Топ-5 по скору:
1. GTCATGGTCATCTGCTACTC | score: 105.0 | GC: 50.0% | strand: +
2. GACAAGTGTGATCACTTGGG | score: 100.0 | GC: 50.0% | strand: +
3. ATGCAGGTGACAGAGACTCT | score: 100.0 | GC: 50.0% | strand: +
4. CAGTTTACACCCGATCCACT | score: 100.0 | GC: 50.0% | strand: +
5. GTAAACTGAGCTTGCTCGCT | score: 100.0 | GC: 50.0% | strand: -
Первый кандидат набрал 105 очков - больше "максимума". Это бонус за гуанин в начале: такие gRNA отлично экспрессируются с U6-промотора.
С флагом --check-offtarget программа проверит каждого кандидата через BLAST. Это займёт несколько минут (NCBI не торопится), но покажет, есть ли в геноме человека опасные похожие участки.
python crispr_finder_v2.py analyze ccr5.fasta --top-n 5 --check-offtarget
============================================================
Ген: CCR5_CDS Homo sapiens C-C motif chemokine receptor 5 (CCR5), mRNA
Длина: 1059 нуклеотидов
============================================================
Найдено кандидатов: 62
Топ-5 по скору:
1. GTCATGGTCATCTGCTACTC | score: 105.0 | GC: 50.0% | strand: +
Проверяю off-target (займёт ~30-60 сек)...
✓ Высокий риск off-target не обнаружен
2. GACAAGTGTGATCACTTGGG | score: 100.0 | GC: 50.0% | strand: +
Проверяю off-target (займёт ~30-60 сек)...
✓ Высокий риск off-target не обнаружен
3. ATGCAGGTGACAGAGACTCT | score: 100.0 | GC: 50.0% | strand: +
Проверяю off-target (займёт ~30-60 сек)...
✓ Высокий риск off-target не обнаружен
4. CAGTTTACACCCGATCCACT | score: 100.0 | GC: 50.0% | strand: +
Проверяю off-target (займёт ~30-60 сек)...
✓ Высокий риск off-target не обнаружен
5. GTAAACTGAGCTTGCTCGCT | score: 100.0 | GC: 50.0% | strand: -
Проверяю off-target (займёт ~30-60 сек)...
✓ Высокий риск off-target не обнаружен
Что можно запилить далее
Мы написали рабочий инструмент, но до уровня Benchling или CRISPOR ему далеко:
ML-скоринг - использовать модели типа Azimuth, обученные на реальных экспериментах
Локальный BLAST - скачать геном и искать без интернета и лимитов NCBI
Позиция мисматчей - мисматч рядом с PAM опаснее, чем на 5'-конце
Визуализация - карта гена с отмеченными сайтами
Batch-режим - много генов за раз с кэшированием
Но для понимания принципов - достаточно.
Итого
За две статьи мы:
Разобрались, как работает CRISPR-Cas9
Написали поиск PAM и извлечение gRNA
Добавили фильтрацию и скоринг
Подключили BLAST для проверки off-target
Протестировали на настоящем гене CCR5
Код в репозитории. Там же FASTA с CCR5.
Послесловие
Как я писал в комментариях к первой части, есть домашние лабораторные наборы с помощью которых можно сделать, например, светящееся пиво или что-то тревожно напоминающее биологическое оружие бактерию выживающую в стрептомицине (антибиотике). Если интересно про что-то такое (и сбоку пайтон) почитать, то голосуйте в опросе под статьёй, устрою в ближайшем будущем (постараюсь).
P.S. Биологам: простите за упрощения. Цель - показать принцип, а не заменить специализированные инструменты.
P.P.S. Не редактируйте людей. Даже если очень хочется. Особенно если очень хочется.
P.P.P.S. Обсудить или покритиковать - велком в телеграм.