Комплексный транскриптомный анализ: от Bulk до Single-cell RNA-Seq

Дизайн биологического эксперимента и технологии подготовки библиотек для RNA-Seq

Геном любой клетки многоклеточного организма практически идентичен и статичен. Транскриптом же напоминает хаотичный, постоянно меняющийся мегаполис: в нейроне и гепатоците активны совершенно разные наборы генов, а их экспрессия меняется каждую минуту в ответ на стресс, температуру или сигналы соседей. Задача RNA-Seq — сделать моментальный снимок этого мегаполиса. Однако то, насколько четким получится этот снимок, зависит не от мощности биоинформатических серверов, а от решений, принятых до того, как первая пробирка коснется льда. Ошибки на этапе пробоподготовки невозможно исправить математическими алгоритмами.

Анатомия биологического эксперимента

Планирование транскриптомного исследования начинается с определения структуры выборки. Главная ошибка на этом этапе — непонимание природы дисперсии (разброса данных) и попытка сэкономить на правильном типе повторностей. Биоинформатические инструменты дифференциальной экспрессии опираются на статистические модели, требующие точной оценки вариативности внутри исследуемых групп.

Биологические и технические повторности

В RNA-Seq мы имеем дело с двумя источниками шума: ошибками самого метода (секвенирования, выделения) и естественной вариативностью живых систем.

Технические повторности возникают, когда берется один и тот же биологический образец (например, гомогенат печени одной конкретной мыши), делится на три пробирки, из них независимо готовятся три библиотеки и секвенируются. Различия между этими тремя результатами покажут техническую погрешность пайплайна. Современные протоколы Illumina, BGI и других платформ настолько точны, что техническая дисперсия минимальна — коэффициенты корреляции между техническими репликами обычно превышают $0.98$ . Делать технические повторности в стандартных RNA-Seq проектах сегодня считается нецелесообразным расходованием бюджета, за исключением случаев валидации совершенно нового кастомного протокола.

Биологические повторности предполагают забор материала от разных независимых организмов (три разные мыши, получавшие одинаковый препарат). Именно они позволяют оценить популяционную дисперсию. Уровень экспрессии гена иммунного ответа TNF у трех генетически идентичных мышей в контрольной группе может отличаться в разы из-за микроскопических различий в их состоянии, микробиоме или перенесенном стрессе. Без знания того, как сильно базово колеблется экспрессия гена у здоровых особей, биоинформатический алгоритм не сможет статистически доказать, что двукратное изменение экспрессии под действием лекарства действительно значимо, а не является случайным выбросом.

Минимальным стандартом для bulk RNA-Seq долгое время считались три биологические повторности на группу ( $N = 3$ ). Однако статистическая мощность такого дизайна крайне низка. Для генов с низкой экспрессией или высокой естественной вариабельностью этого критически мало. Современные рекомендации (например, консорциума ENCODE) требуют $N \geq 4$ для клеточных линий и инбредных животных, и $N \geq 6-10$ для клинических образцов человека, где генетическая гетерогенность пациентов добавляет гигантский слой шума.

Эффект партии (Batch Effect)

Даже при идеальном количестве повторностей эксперимент можно уничтожить неправильным распределением образцов во времени и пространстве. Эффект партии возникает, когда небиологические факторы систематически совпадают с биологическими группами.

Классический сценарий провала: лаборатория исследует опухолевые и здоровые ткани. В понедельник исследователь выделяет РНК из всех контрольных образцов, используя старый набор реактивов. В пятницу он выделяет РНК из опухолей, открыв новую коробку с набором. Через месяц, после секвенирования, анализ главных компонент (PCA) покажет идеальное разделение групп. Эти переменные оказались «сцеплены» (confounded). Математически невозможно определить, вызвана ли разница в профилях экспрессии биологией рака или разницей между наборами реактивов и днями недели.

Правильный дизайн требует рандомизации и блокирования. Если в эксперименте 12 образцов (6 контроль, 6 опыт) и выделение РНК занимает два дня, необходимо в первый день выделить 3 контроля и 3 опыта, и во второй день — оставшиеся 3 контроля и 3 опыта. Аналогично при запуске на секвенаторе: образцы из разных групп должны быть равномерно распределены по дорожкам (lanes) проточной ячейки.

Качество РНК: фундамент библиотеки

Секвенирование РНК начинается с ее выделения, и здесь возникает главная биохимическая проблема: РНК — крайне нестабильная молекула. Одинарная спираль и наличие 2'-OH группы в рибозе делают ее химически уязвимой, а вездесущие ферменты РНКазы способны разрушить образец за минуты.

Для оценки целостности РНК используется капиллярный электрофорез. В тотальной РНК эукариотической клетки более 80% массы приходится на рибосомальную РНК (рРНК), в частности на субъединицы 28S и 18S. В идеальном неповрежденном образце на фореграмме видны два острых высоких пика, причем площадь пика 28S должна быть примерно в два раза больше площади 18S. По мере деградации РНК длинные молекулы рРНК рвутся, пики сглаживаются, и сигнал смещается в сторону коротких фрагментов, образуя «горб» в левой части графика.

Соотношение площадей этих пиков и общий профиль деградации ложатся в основу метрики RIN (RNA Integrity Number). Шкала RIN варьируется от $1.0$ (полностью деградированная РНК) до $10.0$ (идеально интактная РНК).

Показатель RIN диктует, какую технологию подготовки библиотеки допустимо использовать. При $RIN \geq 7$ доступны любые методы. Если RIN падает ниже 5 — что типично для клинических образцов, залитых в парафин (FFPE, Formalin-Fixed Paraffin-Embedded) — стандартные методы приведут к катастрофическим искажениям данных, так как молекулы РНК в таких блоках фрагментированы до кусочков длиной 100-200 нуклеотидов.

Стратегии обогащения: избавление от балласта

Матричная РНК (мРНК), кодирующая белки и представляющая главный интерес для большинства исследователей, составляет всего 1–3% от всей РНК в клетке. Если просто фрагментировать тотальную РНК и отсеквенировать ее, до 95% прочтений (ридов) будут картироваться на гены рибосомальной РНК. Перед созданием библиотеки необходимо провести процедуру обогащения.

Poly-A селекция (Обогащение мРНК)

Большинство зрелых эукариотических мРНК имеют на 3'-конце полиадениновый хвост (Poly-A). Метод использует магнитные шарики, покрытые короткими цепочками олиго-дТ (Oligo-dT). При смешивании шариков с тотальной РНК, поли-А хвосты мРНК гибридизуются с олиго-дТ. Магнит притягивает шарики на стенку пробирки, а вся рибосомальная и транспортная РНК смывается.

Преимущества:

Относительно низкая стоимость реагентов.
Высокая специфичность к белок-кодирующим транскриптам.
Требует меньшей глубины секвенирования (достаточно 20–30 млн ридов на образец), так как секвенируется только информативная фракция.

Уязвимости и ограничения:

Метод абсолютно не работает для прокариот (у бактерий нет стабильных Poly-A хвостов).
Теряются транскрипты, не имеющие Poly-A хвоста. Классический пример — матричные РНК гистоновых белков, которые у большинства эукариот не полиаденилируются. Также теряются многие длинные некодирующие РНК (lncRNA) и кольцевые РНК (circRNA).
3'-bias (смещение к 3'-концу): Это критический нюанс при работе с деградированной РНК. Если длинная молекула мРНК (например, транскрипт гена MUC16 длиной более 10 000 нуклеотидов) разорвана на части из-за низкого RIN, магнитный шарик вытянет только тот фрагмент, на котором остался Poly-A хвост (3'-конец). 5'-конец гена будет безвозвратно утерян. На этапе биоинформатического анализа это проявится как искусственно завышенная экспрессия 3'-участков генов и полное отсутствие прочтений на 5'-концах. Именно поэтому Poly-A селекция требует $RIN \geq 7$ .

Ribo-Zero (Истощение рРНК)

Вместо того чтобы «вытягивать» нужную мРНК, этот метод «выбрасывает» ненужную рРНК (rRNA depletion). В образец добавляются специфические биотинилированные ДНК-зонды, комплементарные последовательностям рибосомальной РНК. Зонды гибридизуются с рРНК, после чего этот комплекс удаляется с помощью стрептавидиновых магнитных шариков (стрептавидин имеет высочайшее сродство к биотину). Все, что осталось в растворе — мРНК, lncRNA, пре-мРНК — идет в библиотеку.

Преимущества:

Позволяет анализировать любые типы РНК, включая некодирующие и гистоновые мРНК.
Спасает образцы с сильной деградацией (FFPE ткани с $RIN < 4$ ). Поскольку метод не полагается на Poly-A хвост, фрагменты мРНК с 5'-конца и из середины гена остаются в растворе и успешно секвенируются, предотвращая 3'-bias.

Уязвимости и ограничения:

В библиотеку попадает много незрелой пре-мРНК с интронами, ядерная РНК и различные виды малых РНК. Это требует увеличения глубины секвенирования (рекомендуется 50–100 млн ридов) для достижения той же статистической мощности при подсчете экспрессии экзонов, так как часть ридов «уйдет» на интронные последовательности.

Сравнительная таблица стратегий

Характеристика	Poly-A селекция	Ribo-Zero (rRNA Depletion)
Целевая фракция	Только РНК с Poly-A хвостом	Вся РНК, кроме рибосомальной
Требование к качеству (RIN)	Строгое ( $RIN \geq 7$ )	Мягкое (подходит для деградированной РНК)
Захват lncRNA и гистонов	Плохой (только полиаденилированные)	Отличный
Наличие интронных ридов	Минимальное	Высокое (захватывает пре-мРНК)
Организмы	Только эукариоты	Эукариоты и прокариоты

Синтез кДНК и проблема направленности (Strandedness)

Секвенаторы Illumina не умеют читать РНК напрямую — им нужна двухцепочечная ДНК. Поэтому после обогащения и химической фрагментации РНК (обычно до кусочков в 200–300 нуклеотидов) необходимо провести обратную транскрипцию.

Исторически первые протоколы RNA-Seq создавали обычную двухцепочечную кДНК (комплементарную ДНК). К ней пришивались адаптеры, и она отправлялась на секвенирование. Проблема заключалась в потере информации о том, с какой именно цепи геномной ДНК был считан транскрипт (с «плюс» или «минус» цепи).

В плотных эукариотических геномах часто встречаются перекрывающиеся гены. Например, ген NR1D1 закодирован на одной цепи ДНК, а ген THRA — на противоположной, при этом их 3'-концы физически перекрываются. Если прочитать фрагмент из зоны перекрытия классическим не-направленным методом, невозможно определить, продуктом какого из двух генов является этот рид.

Современным стандартом является направленный (stranded) RNA-Seq, чаще всего реализуемый через метод включения dUTP.

Механика метода элегантна в своей биохимической простоте:

Синтез первой цепи: Обратная транскриптаза строит первую цепь кДНК на матрице РНК, используя стандартные дезоксирибонуклеотиды (dATP, dCTP, dGTP, dTTP).
Синтез второй цепи: РНК-матрица удаляется ферментом РНКазой H. ДНК-полимераза строит вторую цепь кДНК. Ключевой трюк заключается в том, что в реакционную смесь вместо тимина (dTTP) добавляют урацил (dUTP). Вторая цепь оказывается «помеченной» урацилами.
Лигирование адаптеров: К двухцепочечным фрагментам пришиваются Y-образные адаптеры секвенирования.
Разрушение метки: Перед финальной ПЦР-амплификацией в пробирку добавляют фермент урацил-ДНК-гликозилазу (UDG/USER), который специфично распознает и вырезает урацил, разрушая вторую цепь.
Амплификация: Вторая цепь уничтожена. ПЦР идет только с первой цепи, которая точно соответствует исходной РНК.

В результате биоинформатик получает риды, которые картируются на геном со строгим указанием исходной цепи. Это радикально повышает точность подсчета экспрессии в сложных локусах и позволяет выявлять антисмысловые РНК (antisense RNA), играющие важную роль в регуляции транскрипции.

Параметры секвенирования: глубина и длина рида

Финальный этап дизайна — заказ параметров на секвенаторе. Необходимо определить глубину (количество прочтений на образец) и режим чтения.

Глубина секвенирования (Sequencing Depth) определяет чувствительность эксперимента. Человеческий геном содержит около 20 000 белок-кодирующих генов. Если цель — найти дифференциально экспрессируемые гены со средним и высоким уровнем транскрипции, 20–30 миллионов ридов на образец (для Poly-A библиотек) обеспечивают достаточное покрытие. Если задача — детектировать редкие транскрипты, факторы транскрипции или анализировать альтернативный сплайсинг, глубину увеличивают до 50–100 миллионов. Дальнейшее увеличение глубины подчиняется закону убывающей отдачи: затраты растут кратно, но открываются лишь единичные новые низкокопийные гены, статистическая значимость которых часто сомнительна.

Режим чтения: Single-end (SE) против Paired-end (PE). В режиме SE секвенатор читает фрагмент кДНК только с одного конца (обычно 50–75 нуклеотидов). В режиме PE прибор читает фрагмент с обоих концов (например, $2 \times 150$ нуклеотидов).

Режим SE дешевле и полностью покрывает базовую задачу — подсчет уровня экспрессии известных генов (Gene Counting). Короткого участка в 50 нуклеотидов достаточно, чтобы алгоритм картирования (например, STAR) однозначно определил, из какого гена выпал этот рид.

Режим PE обязателен в следующих случаях:

Изучение альтернативного сплайсинга (риды с двух концов длинного фрагмента могут попасть в разные экзоны, разделенные длинным интроном, доказывая их фактическое соединение в зрелой мРНК).
Сборка транскриптома de novo (без референсного генома), где парные концы помогают алгоритмам (например, Trinity) строить длинные контиги.
Работа с аллель-специфичной экспрессией, где требуется максимальное покрытие полиморфизмов (SNP) на одной физической молекуле.

Первичная обработка сырых данных: контроль качества и путь от FASTQ к матрице экспрессии

После нескольких недель планирования эксперимента, выделения РНК и ожидания очереди на секвенатор, лаборатория получает ссылку на скачивание. По ссылке — архив весом в десятки гигабайт, содержащий сотни миллионов коротких текстовых строк из четырех букв: A, T, G, C. На этом этапе биология временно отступает на второй план, уступая место вычислительной математике и алгоритмам обработки текстов. Задача первичного анализа — превратить этот массив разрозненных фрагментов генома в компактную математическую матрицу, где для каждого гена будет указано точное количество его копий в каждом образце.

Анатомия сырых данных: формат FASTQ

Секвенаторы платформ Illumina выдают результаты в текстовом формате FASTQ. В отличие от классического формата FASTA, который хранит только саму нуклеотидную последовательность, FASTQ содержит критически важный дополнительный слой информации — уверенность машины в правильности прочтения каждой конкретной буквы.

Каждый прочитанный фрагмент кДНК (рид) представлен в файле ровно четырьмя строками.

@ERR1234567.1 1 length=100
GATCTGATAACTCGGTCGAAATTTTCAAGCGTTATAG...
+
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII...

Первая строка всегда начинается с символа @ и содержит уникальный идентификатор рида, координаты кластера на проточной ячейке секвенатора и техническую информацию о запуске. Вторая строка — это непосредственно прочитанные нуклеотиды (A, C, G, T и символ N, если оптическая система не смогла распознать сигнал). Третья строка выступает разделителем и содержит символ +. Четвертая строка — строка качества, длина которой строго равна длине нуклеотидной последовательности во второй строке.

Для оценки качества используется логарифмическая метрика Phred score, обозначаемая как $Q$ . Она связывает вероятность ошибочного распознавания нуклеотида, обозначаемую как $P$ , с целым числом:

Q = -10 \log_{10} P

В этой формуле $P$ представляет собой вероятность того, что секвенатор ошибся при вызове базы (base calling). Если вероятность ошибки составляет $1$ к $1000$ ( $P = 0.001$ ), то показатель качества $Q$ будет равен $30$ .

Значение $Q$	Вероятность ошибки ( $P$ )	Точность прочтения
10	1 к 10 ( $0.1$ )	90.0%
20	1 к 100 ( $0.01$ )	99.0%
30	1 к 1000 ( $0.001$ )	99.9%
40	1 к 10000 ( $0.0001$ )	99.99%

В биоинформатике негласным стандартом приемлемого качества считается $Q \geq 20$ , а отличного — $Q \geq 30$ . Чтобы уместить двузначные числа Phred score в один символ текстового файла и не нарушить структуру, используется кодировка ASCII (стандарт Phred+33). К значению $Q$ прибавляется $33$ , и полученное число интерпретируется как код символа в таблице ASCII. При $Q=40$ в файл запишется символ с кодом $73$ , что соответствует заглавной латинской букве I.

Контроль качества и «стрижка» ридов (Trimming)

Доверие к сырым данным без предварительной проверки часто приводит к ложноположительным результатам при поиске дифференциально экспрессирующихся генов. Первый шаг — агрегация метрик качества миллионов ридов с помощью программ оценки качества (FastQC) и их объединение в единый отчет (MultiQC).

Две главные проблемы, выявляемые на этом этапе:

Падение качества на 3'-конце рида. Химия секвенирования путем синтеза (Sequencing by Synthesis, SBS) подвержена постепенному накоплению ошибок. Процесс опирается на синхронное присоединение флуоресцентно-меченых нуклеотидов тысячами идентичных молекул в одном кластере. С каждым новым циклом часть молекул отстает (phasing), а часть забегает вперед (pre-phasing). Флуоресцентный сигнал становится всё более «размытым». В результате первые 50–70 нуклеотидов могут иметь средний Phred score около $35$ , а последние 10–20 нуклеотидов проваливаются ниже $20$ .
Контаминация адаптерами. Библиотеки для секвенирования содержат синтетические адаптеры на концах фрагментов. Если исходный биологический фрагмент кДНК оказался короче, чем заданная длина чтения (например, фрагмент 100 пар оснований при режиме секвенирования 150 пар оснований), полимераза прочитает весь биологический фрагмент и продолжит синтез, читая адаптер на противоположном конце (read-through). Попытка найти эту химерную последовательность в референсном геноме обречена на провал.

Для решения этих проблем применяется процедура «стрижки» (trimming) с помощью инструментов типа Trimmomatic или Cutadapt. Программа сканирует каждый рид и выполняет две операции. Сначала она ищет точные совпадения с известными последовательностями адаптеров Illumina и отрезает их. Затем применяется алгоритм скользящего окна (sliding window) для удаления низкокачественных концов.

Алгоритм скользящего окна анализирует не каждый нуклеотид по отдельности, а их группы. Окно заданного размера движется от 3'-конца к 5'-концу. Если среднее качество нуклеотидов в окне падает ниже порогового значения, рид обрезается в этой точке.

Например, задано окно размером 4 нуклеотида и порог $Q=20$ . Программа берет четыре нуклеотида с 3'-конца со значениями Phred [22, 18, 15, 10]. Среднее значение равно $16.25$ . Поскольку $16.25 < 20$ , эти нуклеотиды отсекаются. Окно сдвигается левее. Если следующие четыре нуклеотида имеют качество [30, 25, 22, 18], их среднее равно $23.75$ . Это значение выше порога, обрезка останавливается, и оставшаяся часть рида сохраняется. Если после всех процедур рид становится слишком коротким (менее 36 пар оснований), он полностью удаляется из датасета.

Сплайсированное выравнивание: поиск координат на геноме

Очищенные текстовые строки необходимо привязать к конкретным генам. Этот процесс называется картированием или выравниванием (alignment). В анализе ДНК задача сводится к поиску наиболее похожего участка в референсном геноме. Однако в транскриптомике эукариот возникает фундаментальная биологическая преграда — сплайсинг.

Зрелая мРНК, из которой строится библиотека для секвенирования, не содержит интронов. Но референсный геном, к которому мы прикладываем риды, содержит их в полном объеме. Если рид случайно попал на границу двух экзонов, его первая половина должна прикрепиться к одному участку генома, а вторая — к другому, находящемуся на расстоянии десятков или сотен тысяч пар нуклеотидов. Стандартные ДНК-выравниватели (например, BWA или Bowtie2) штрафуют алгоритм за наличие длинных пробелов (гэпов). Они воспримут такой рид как содержащий гигантскую делецию и просто отбросят его как мусорный.

Для RNA-Seq требуются специализированные «сплайсинг-ориентированные» алгоритмы, самым популярным из которых является STAR (Spliced Transcripts Alignment to a Reference).

Алгоритм STAR использует стратегию Maximum Mappable Prefix (MMP). Он берет начало рида и ищет в геноме самую длинную последовательность, которая совпадает с ним идеально. Как только алгоритм натыкается на несовпадение (границу экзона), поиск MMP останавливается. STAR фиксирует эту координату и начинает искать идеальное совпадение для оставшейся части рида (суффикса) в других участках генома. Если вторая часть рида находится ниже по течению, а на границах разрыва в геноме обнаруживаются канонические мотивы сплайсинга (GT на 5'-конце интрона и AG на 3'-конце), алгоритм делает вывод: рид пересекает интрон.

Форматы SAM/BAM и язык CIGAR

Результат работы выравнивателя сохраняется в формате SAM (Sequence Alignment Map) или его бинарной, сжатой версии BAM. Это огромная таблица, где для каждого исходного рида указана хромосома, точная координата начала выравнивания, качество картирования (MAPQ) и специфическая CIGAR-строка.

CIGAR (Compact Idiosyncratic Gapped Alignment Report) — это буквенно-цифровой код, описывающий геометрию выравнивания рида относительно референса.

Ключевые операторы CIGAR:

M (Alignment Match): нуклеотид рида сопоставлен с нуклеотидом генома. Оператор M не гарантирует биохимическую идентичность букв, он лишь указывает на геометрическое выравнивание в этой позиции (там может находиться однонуклеотидный полиморфизм, SNP).
I (Insertion): в риде есть нуклеотиды, которых нет в референсе (вставка).
D (Deletion): в референсе есть нуклеотиды, которых нет в риде (делеция).
N (Skipped region): специфический для RNA-Seq оператор, обозначающий интрон.
S (Soft clipping): нуклеотиды присутствуют в риде, но алгоритм решил их проигнорировать при выравнивании (часто возникает, если на конце рида остался кусок адаптера, который не удалось удалить на этапе тримминга).

Рассмотрим три сценария для рида длиной 100 нуклеотидов:

Рид идеально лег на один экзон без мутаций и разрывов. Его CIGAR: 100M.
Рид пересек границу сплайсинга: 40 нуклеотидов попали в первый экзон, затем следует интрон длиной 5000 баз, а оставшиеся 60 нуклеотидов легли на второй экзон. CIGAR: 40M5000N60M. Именно оператор N позволяет биоинформатикам реконструировать альтернативный сплайсинг.
Первые 90 нуклеотидов совпали с геномом, а последние 10 оказались остатком адаптера. Алгоритм «отрежет» их виртуально. CIGAR: 90M10S.

Альтернативный путь: псевдовыравнивание

Сборка BAM-файлов с помощью STAR требует значительных вычислительных мощностей (от 30 ГБ оперативной памяти для генома человека). Если цель исследования — только оценка уровня экспрессии известных генов, а не поиск новых изоформ или мутаций, применяется альтернативный подход: псевдовыравнивание (инструменты Salmon, Kallisto).

Псевдовыравниватели картируют риды не на полный геном с интронами, а на транскриптом — базу данных уже известных последовательностей зрелых мРНК. Вместо посимвольного выравнивания они разбивают риды и референсные транскрипты на короткие фрагменты фиксированной длины, называемые k-мерами.

Если последовательность рида — ATCGTA, а размер $k$ равен 3, то алгоритм разобьет его на четыре 3-мера: ATC, TCG, CGT, GTA. Из этих k-меров строится граф де Брюйна. Алгоритм не пытается выяснить точную геномную координату рида. Он лишь отвечает на вопрос: «С каким известным транскриптом этот набор k-меров совместим?». Отказ от точного позиционирования снижает требования к памяти до 4–8 ГБ и ускоряет процесс в десятки раз, выдавая на выходе сразу готовые оценки экспрессии без промежуточных BAM-файлов.

Квантификация: от координат к матрице экспрессии

Если использовался классический путь через геномное выравнивание (STAR), финальным этапом первичной обработки становится квантификация — подсчет количества ридов, попавших в границы каждого гена. Стандартом индустрии для этой задачи является программа featureCounts (из пакета Subread).

На вход программе подаются BAM-файлы с координатами ридов и аннотация генома (GTF/GFF файл) — справочник, в котором указано, на какой хромосоме и в каких координатах начинаются и заканчиваются экзоны всех известных генов.

Процесс подсчета сталкивается с несколькими пространственными коллизиями, требующими жестких правил разрешения:

Перекрытие генов на разных цепях. На участке хромосомы может находиться ген A на прямой цепи (+) и ген B на обратной цепи (-). Если рид попал в зону их перекрытия, а библиотека готовилась по ненаправленному протоколу, определить источник рида невозможно. featureCounts по умолчанию отбросит его (присвоит статус ambiguous). Если же применялся направленный dUTP-протокол, параметр strandedness укажет программе засчитать рид только тому гену, с чьей цепи он был транскрибирован.
Мульти-маппинг (Multi-mapping). Некоторые риды идеально выравниваются в нескольких местах генома одновременно. Это типично для семейств паралогичных генов или псевдогенов. В стандартном bulk RNA-Seq анализе такие риды принято исключать из подсчета, чтобы не создавать искусственное завышение экспрессии в повторяющихся участках.
Парные чтения (Paired-end). Если секвенирование проводилось с двух концов фрагмента, featureCounts настраивается на подсчет не отдельных ридов, а целых фрагментов (параметр countReadPairs). Если левый и правый риды корректно легли на один ген, это засчитывается как один фрагмент, а не как два независимых события.

Результатом работы featureCounts является таблица — матрица сырых каунтов (raw count matrix). В ней строки соответствуют генам, столбцы — биологическим образцам (контроль 1, контроль 2, опыт 1...), а на пересечении стоят целые числа: сколько уникальных фрагментов кДНК было приписано данному гену в данном образце.

Эта матрица является водоразделом в транскриптомном анализе. Физические молекулы, химические реакции секвенирования, оптические сигналы, текстовые строки и геометрические координаты окончательно свернулись в абстрактную математическую форму. Дальнейшая работа будет опираться исключительно на статистические распределения этих чисел для поиска биологического смысла.

Статистический анализ bulk RNA-Seq: методы нормализации и оценка дифференциальной экспрессии генов

В матрице каунтов, полученной после выравнивания, ген ACTB имеет 500 прочтений в контрольном образце и 1000 в опытном. Означает ли это, что экспрессия гена выросла ровно в два раза? Нет. Возможно, секвенатор сгенерировал в два раза больше данных для опытного образца. А возможно, в опытном образце «замолчал» другой ген, который в норме забирал на себя половину всех прочтений, и ACTB просто занял освободившееся место на проточной ячейке, хотя его реальное количество в клетке не изменилось. Прямое сравнение сырых прочтений (raw counts) между образцами математически бессмысленно и ведет к ложным биологическим выводам, которые впоследствии могут направить целое исследование по ложному пути.

Две проблемы сырых каунтов: глубина и композиция

Чтобы сделать данные сопоставимыми, их необходимо нормализовать. В экспериментах RNA-Seq мы сталкиваемся с двумя фундаментальными искажениями, которые нужно компенсировать до начала любого статистического тестирования.

Первое искажение — глубина секвенирования (sequencing depth). Разные библиотеки почти всегда получают разное количество тотальных прочтений. Это связано с микроскопическими различиями в концентрации ДНК при загрузке на проточную ячейку (flow cell) секвенатора. Если библиотека А содержит 20 миллионов ридов, а библиотека Б — 40 миллионов, то при прочих равных условиях любой ген в библиотеке Б получит в два раза больше каунтов. Это технический артефакт, не имеющий отношения к биологии.

Второе, более коварное искажение — эффект композиции РНК (RNA composition bias). Проточная ячейка секвенатора имеет фиксированную емкость. Секвенирование — это процесс случайной выборки молекул из пула (sampling).

Рассмотрим упрощенную модель клетки, в которой экспрессируются всего три гена. В контрольном состоянии ген А производит 100 транскриптов, ген Б — 100, ген В — 10. Всего в клетке 210 молекул РНК. Мы вводим токсин, который вызывает колоссальную экспрессию гена В, так что он начинает производить 1000 транскриптов. Гены А и Б никак не отреагировали на токсин, их абсолютное количество осталось прежним (по 100 транскриптов). Теперь в клетке 1200 молекул РНК.

Если мы секвенируем оба образца с одинаковой глубиной, скажем, по 210 ридов на каждый, мы получим следующую картину:

В контроле: Ген А получит ~100 ридов, Ген Б ~100 ридов.
В опыте: Ген В заберет на себя подавляющее большинство ридов (около $1000 / 1200 \approx 83\%$ ). На гены А и Б останется суммарно около 34 ридов (по 17 на каждый).

Если мы просто разделим каунты каждого гена на общее число ридов в библиотеке, нам покажется, что экспрессия генов А и Б резко упала (со 100 до 17). На самом деле их абсолютное количество в клетке не изменилось, они просто были вытеснены из выборки одним супер-экспрессирующимся транскриптом.

Почему TPM и RPKM не подходят для дифференциальной экспрессии

Исторически первыми методами нормализации были метрики RPKM (Reads Per Kilobase of transcript per Million mapped reads) и TPM (Transcripts Per Million). Они пытались решить сразу две задачи: учесть длину гена (длинные гены при случайной фрагментации генерируют больше фрагментов, чем короткие) и учесть глубину библиотеки.

Сегодня стандартом для оценки относительной представленности транскрипта внутри одного образца считается TPM. Формула расчета TPM выглядит так:

TPM_i = \left( \frac{q_i / l_i}{\sum_j (q_j / l_j)} \right) \times 10^6

Где:

$TPM_i$ — нормализованное значение для гена $i$ .
$q_i$ — количество сырых прочтений (counts), картированных на ген $i$ .
$l_i$ — эффективная длина гена $i$ в килобазах.
$\sum_j (q_j / l_j)$ — сумма отношений каунтов к длине для всех генов $j$ в данном образце.

Суть TPM в том, что сумма всех значений TPM в любом образце всегда равна ровно $10^6$ (одному миллиону). Это позволяет корректно сравнивать гены внутри одного образца: если TPM гена А равен 100, а гена Б равен 50, значит, в пуле мРНК физических молекул гена А в два раза больше. Эта метрика будет крайне полезна нам в будущем, при переходе к single-cell RNA-Seq, где TPM-подобные трансформации часто используются для кластеризации клеток.

Однако для поиска дифференциально экспрессирующихся генов (DGE) между разными образцами TPM категорически не подходит. Из-за того, что сумма всегда равна миллиону, метрика TPM жестко подвержена описанному выше эффекту композиции. Если один ген заберет на себя 800 000 TPM, остальные 20 000 генов будут вынуждены делить между собой оставшиеся 200 000, искусственно занижая свои показатели. Сравнение TPM между контролем и опытом приведет к массовому обнаружению ложно-подавленных генов.

Алгоритм Median of Ratios (DESeq2)

Для корректного сравнения образцов друг с другом требуется метод, устойчивый к выбросам (robust). Стандартом индустрии стал алгоритм Median of Ratios, реализованный в пакете DESeq2. Его главная идея — найти масштабирующий коэффициент (size factor) для каждой библиотеки, опираясь на гены, экспрессия которых стабильна во всех условиях.

Алгоритм работает в четыре шага:

Создание псевдо-референсного образца. Для каждого гена вычисляется среднее геометрическое его каунтов по всем образцам эксперимента. Среднее геометрическое используется потому, что оно менее чувствительно к экстремальным значениям, чем среднее арифметическое. Если в каком-то образце ген имеет 0 каунтов, среднее геометрическое обнуляется, и ген исключается из расчетов size factor (чтобы избежать деления на ноль в следующем шаге).
Расчет отношений. Для каждого гена в каждом образце его сырой каунт делится на значение из псевдо-референса. Мы получаем массив отношений.
Вычисление медианы. Для каждого образца берется медиана всех полученных отношений. Эта медиана и есть size factor (масштабирующий множитель) данной библиотеки.
Нормализация. Сырые каунты каждого образца делятся на его size factor.

Почему используется именно медиана, а не среднее арифметическое? Медиана абсолютно нечувствительна к экстремальным выбросам. Возвращаясь к нашему примеру с токсином: если ген В вырос в 10 раз, он сильно исказит среднее арифметическое отношений. Но медиана массива отношений останется на уровне тех генов (housekeeping genes), чья экспрессия не изменилась (их отношение к псевдо-референсу будет близко к 1). Таким образом, Median of Ratios элегантно решает проблему эффекта композиции, автоматически игнорируя гены с аномальными изменениями при расчете поправки на глубину секвенирования.

Статистическая природа данных: Отрицательное биномиальное распределение

Получив нормализованные каунты, мы не можем просто применить классический t-критерий Стьюдента для оценки статистической значимости различий.

Во-первых, каунты РНК-секвенирования — это дискретные величины (целые числа: 0, 1, 2, 3...), а не непрерывные, для которых создано нормальное распределение. Во-вторых, в RNA-Seq дисперсия (разброс значений) зависит от среднего значения. У низкоэкспрессируемых генов разброс между повторностями небольшой в абсолютных цифрах (например, 5, 7 и 4 каунта), а у высокоэкспрессируемых — огромный (например, 10000, 12000 и 9000 каунтов).

Изначально в ранних методах предполагалось, что каунты подчиняются распределению Пуассона, где дисперсия строго равна среднему. Если мы секвенируем одну и ту же библиотеку несколько раз (технические повторности), данные действительно подчиняются Пуассону — разброс обусловлен только случайностью попадания молекулы в пору секвенатора.

Однако биологические системы обладают сверхдисперсией (overdispersion). Различия между биологическими повторностями (разными мышами, разными клеточными культурами, разными пациентами) всегда превышают техническую погрешность секвенатора. У пациентов разный возраст, генетический фон, уровень стресса — все это вносит огромный вклад в вариативность экспрессии.

Поэтому DESeq2 и аналогичные пакеты (например, edgeR) моделируют данные с помощью отрицательного биномиального распределения (Negative Binomial distribution, NB). Связь между дисперсией и средним в NB-распределении описывается формулой:

Var(K) = \mu + \alpha \mu^2

Где:

$Var(K)$ — ожидаемая дисперсия каунтов для гена.
$\mu$ — среднее значение экспрессии гена.
$\alpha$ — параметр дисперсии (dispersion parameter), отражающий биологическую вариабельность.

Если $\alpha = 0$ , квадратичный член исчезает, формула превращается в $Var(K) = \mu$ , что в точности соответствует распределению Пуассона (только технический шум). Но в реальных биологических данных $\alpha > 0$ , и квадратичный член $\alpha \mu^2$ позволяет корректно моделировать сильный разброс у высокоэкспрессируемых транскриптов, учитывая индивидуальные различия между живыми организмами.

Проблема малых выборок и сжатие дисперсии (Shrinkage)

Чтобы статистический тест (в DESeq2 используется тест Вальда) определил, значимо ли изменение экспрессии, ему необходимо точно знать дисперсию гена. Но в типичном RNA-Seq эксперименте у нас крайне мало данных — часто всего 3 или 4 биологические повторности на группу. Оценить параметр $\alpha$ по трем точкам с приемлемой точностью математически невозможно — оценка будет крайне нестабильной. У одних генов дисперсия случайно окажется заниженной (и мы получим ложноположительный результат, так как алгоритм решит, что ген очень стабилен), у других — завышенной.

DESeq2 решает эту проблему с помощью техники сжатия дисперсии (dispersion shrinkage), опираясь на парадигму Эмпирического Байеса и заимствования информации (borrowing information) между генами.

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

Сначала рассчитывается индивидуальная оценка дисперсии для каждого из 20 000 генов по отдельности.
Затем через это облако точек проводится регрессионная кривая — она показывает ожидаемую дисперсию для любого заданного уровня экспрессии.
Наконец, индивидуальные оценки дисперсии "подтягиваются" (shrink) к этой кривой.

Если ген имел аномально низкую дисперсию исключительно из-за случайности малой выборки, алгоритм принудительно увеличит ее до уровня кривой. Это защищает исследователя от ложных открытий (false positives), делая статистический тест более консервативным и надежным.

Исключение делается только для генов с экстремально высокой истинной дисперсией (например, генов иммунного ответа, которые могут сильно варьировать даже в норме). Если индивидуальная дисперсия гена находится слишком высоко над кривой, алгоритм распознает его как истинный биологический выброс и не сжимает его значение вниз, чтобы не спровоцировать ложноположительный результат при сравнении групп.

Множественное тестирование и контроль FDR

Когда модель построена и финальные дисперсии оценены, пакет вычисляет $p$ -value для каждого гена, проверяя нулевую гипотезу: «изменение экспрессии этого гена между контрольной и опытной группами равно нулю».

Здесь возникает фундаментальная проблема анализа омиксных данных — проблема множественного тестирования. В геноме человека около 20 000 белок-кодирующих генов. Если мы установим стандартный порог значимости $\alpha = 0.05$ (вероятность ошибки I рода 5%), то даже при сравнении двух абсолютно идентичных образцов (где нет никакой биологической разницы), статистика выдаст нам $20 000 \times 0.05 = 1000$ генов со значимым $p$ -value исключительно за счет случайных флуктуаций.

Чтобы не утонуть в ложноположительных результатах, сырые $p$ -value необходимо скорректировать. Жесткая поправка Бонферрони (умножение $p$ -value на число тестов) здесь не подходит — она уничтожит все реальные биологические сигналы, так как тестов слишком много. В транскриптомике золотым стандартом является контроль FDR (False Discovery Rate) с помощью поправки Бенджамини-Хохберга.

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

Математика поправки Бенджамини-Хохберга элегантна и опирается на ранжирование. Гены сортируются по возрастанию их сырого $p$ -value. Затем для каждого гена вычисляется скорректированное значение ( $p_{adj}$ ) по формуле:

p_{adj} = p \times \frac{m}{k}

Где:

$p_{adj}$ — скорректированное $p$ -value (оценка FDR).
$p$ — исходное $p$ -value гена.
$m$ — общее количество протестированных генов.
$k$ — ранг гена в отсортированном списке (от 1 — самый значимый, до $m$ — наименее значимый).

Ген	Сырое $p$ -value	Ранг ( $k$ )	Всего генов ( $m$ )	Расчет $p_{adj}$	Итог $p_{adj}$
Ген X	0.0001	1	10000	$0.0001 \times (10000 / 1)$	1.0000 (ограничивается сверху) -> 0.0001
Ген Y	0.0010	2	10000	$0.0010 \times (10000 / 2)$	0.0050
Ген Z	0.0500	100	10000	$0.0500 \times (10000 / 100)$	5.0000 (ограничивается 1.0)

Примечание: в реальных алгоритмах $p_{adj}$ дополнительно корректируется, чтобы значения не убывали при движении вниз по списку, и никогда не превышали 1.0.

Если мы отфильтруем результаты по $p_{adj} < 0.05$ , это будет означать: «среди полученного списка дифференциально экспрессирующихся генов мы готовы смириться с тем, что 5% являются статистическим шумом, но 95% — реальные биологические изменения».

Log2 Fold Change и независимая фильтрация

Статистическая значимость ( $p_{adj}$ ) говорит нам лишь о том, насколько мы уверены, что изменение отлично от нуля. Но она ничего не говорит о масштабе этого изменения. Ген с гигантской базовой экспрессией и нулевой дисперсией может иметь $p_{adj} = 10^{-10}$ , но его экспрессия изменилась всего на 5%. Биологический смысл такого изменения часто равен нулю — клетка не заметит разницы.

Для оценки масштаба используется Log2 Fold Change (LFC) — логарифм по основанию 2 от отношения нормализованной экспрессии в опыте к контролю. Почему именно логарифм по основанию 2? Он делает изменения симметричными. Если ген вырос в 2 раза, отношение равно 2, а $LFC = 1$ . Если ген упал в 2 раза, отношение равно 0.5, а $LFC = -1$ . Без логарифмирования падение экспрессии сжималось бы в узком диапазоне от 0 до 1, а рост уходил бы в бесконечность, что делает визуализацию (например, на Volcano plot) невозможной.

На практике исследователи применяют двойной фильтр: ищут гены, у которых одновременно $p_{adj} < 0.05$ и $|LFC| > 1$ (изменение более чем в два раза в любую сторону).

Важным этапом оптимизации статистической мощности перед применением поправки Бенджамини-Хохберга является независимая фильтрация (independent filtering). Гены с очень низким числом прочтений (например, в среднем менее 10 каунтов на образец) обладают настолько высоким уровнем пуассоновского шума, что алгоритм математически не сможет признать их изменения статистически значимыми, даже если LFC велик.

Оставляя такие шумовые гены в анализе, мы лишь увеличиваем параметр $m$ (общее число тестов) в формуле Бенджамини-Хохберга. Чем больше $m$ , тем сильнее штраф для всех остальных генов. Пакет DESeq2 автоматически отсекает такие малоинформативные гены до этапа множественного тестирования. Это уменьшает $m$ и делает поправку менее жесткой для хороших, высокоэкспрессируемых генов, позволяя найти больше реальных биологических сигналов.

Функциональная аннотация и интерпретация результатов: анализ обогащения GO, KEGG и GSEA

Вы смотрите на итоговую таблицу после запуска DESeq2 и применения поправок на множественное тестирование. В ней 1500 строк — дифференциально экспрессирующихся генов (DEGs). Смотреть на этот список глазами и пытаться уловить биологический смысл — занятие бессмысленное. Даже если в топе таблицы вы узнаете знакомые гены-регуляторы вроде TP53 или TNF, человеческий мозг не способен синтезировать из тысяч разрозненных строк целостную картину клеточных изменений. Клетка не меняет экспрессию случайного набора генов; она активирует или подавляет скоординированные функциональные модули: метаболические пути, каскады передачи сигнала, структурные белковые комплексы.

Переход от списка индивидуальных генов к пониманию этих глобальных модулей называется функциональной аннотацией и анализом обогащения (Enrichment Analysis). Это мост между сухой статистикой матриц каунтов и реальной биологией, позволяющий делать обоснованные выводы о том, почему клетка ведет себя именно так в ответ на стимул или болезнь.

Базы данных: формализация биологического знания

Чтобы биоинформатический алгоритм мог найти закономерности в списке генов, биологическое знание должно быть машиночитаемым. Десятилетиями ученые читали статьи и вручную собирали информацию о функциях генов, формируя специализированные базы данных. Эти базы группируют гены в осмысленные наборы (Gene Sets). Для успешного анализа необходимо понимать архитектуру трех фундаментальных систем: Gene Ontology, KEGG и MSigDB.

Gene Ontology (GO): строгая иерархия терминов

Gene Ontology — это не просто плоский список категорий. Это строгий, универсальный для всех организмов словарь, описывающий свойства генов. GO разделена на три независимые ветви (домены), каждая из которых отвечает на свой вопрос:

Biological Process (BP) — биологический процесс. Отвечает на вопрос «в какой глобальной задаче участвует продукт гена?». Примеры: трансляция, клеточное деление, апоптоз, иммунный ответ.
Molecular Function (MF) — молекулярная функция. Отвечает на вопрос «что физически делает молекула на биохимическом уровне?». Примеры: киназная активность, связывание АТФ, ДНК-хеликазная активность.
Cellular Component (CC) — клеточный компонент. Отвечает на вопрос «где локализован продукт гена в клетке?». Примеры: митохондриальная внутренняя мембрана, рибосома, ядро, синапс.

Главная математическая особенность GO заключается в ее архитектуре. Это направленный ациклический граф (Directed Acyclic Graph, DAG). Термины организованы от широких понятий к узким, при этом один термин может иметь несколько «родителей».

Правило истинного пути (True Path Rule): если ген аннотирован узким термином, он автоматически считается аннотированным и всеми родительскими терминами вплоть до корня графа.

Gene Ontology Consortium

Если ген BRCA1 аннотирован термином «репарация двунитевых разрывов ДНК путем гомологичной рекомбинации», алгоритм автоматически припишет этот ген к терминам «репарация ДНК», «реакция на стресс» и базовому «биологический процесс». Эта избыточность обеспечивает невероятную полноту данных, но приводит к тому, что в результатах анализа часто появляются десятки математически значимых, но биологически дублирующих друг друга терминов.

KEGG и Reactome: метаболические и сигнальные пути

В отличие от GO, которая классифицирует гены по изолированным свойствам, базы данных путей (Pathway Databases) фокусируются на взаимодействиях.

База KEGG (Kyoto Encyclopedia of Genes and Genomes) представляет собой собранные вручную карты, показывающие, как молекулы взаимодействуют друг с другом в рамках конкретного пути (например, «Гликолиз» или «Сигнальный путь Wnt»). Если GO говорит нам, что гены A, B и C участвуют в апоптозе, то KEGG показывает топологию: белок A фосфорилирует белок B, который перемещается в ядро и ингибирует транскрипцию гена C.

Альтернативой KEGG выступает база Reactome. Она отличается более высокой детализацией биохимических реакций и чаще обновляется сообществом. В классическом анализе обогащения топология путей игнорируется (путь рассматривается просто как мешок генов), но для продвинутых методов (таких как Signaling Pathway Impact Analysis, SPIA) направления связей между узлами графа критически важны для оценки того, активирован путь или подавлен.

MSigDB: золотой стандарт для GSEA

Molecular Signatures Database (MSigDB) — это курируемая коллекция наборов генов, изначально созданная для метода GSEA. Самая важная ее часть — коллекция Hallmark (H). Она содержит 50 специфических наборов генов, которые описывают фундаментальные биологические состояния (например, гипоксия, ответ на интерферон гамма, переход от G2 к M фазе клеточного цикла). Наборы Hallmark очищены от избыточности и шума, что делает их идеальной стартовой точкой для интерпретации транскриптомов млекопитающих.

Over-Representation Analysis (ORA): поиск аномальных пересечений

Исторически первым, самым быстрым и интуитивным методом интерпретации стал Over-Representation Analysis (ORA). Его логика опирается на пересечение множеств: мы берем список генов, которые достоверно изменили экспрессию (например, прошли жесткие фильтры $LFC > 1$ и $FDR < 0.05$ ), и проверяем, нет ли среди них аномально высокой доли генов из какого-то известного биологического пути.

Математически ORA использует гипергеометрическое распределение или точный тест Фишера. Классическая аналогия для этого теста — задача про урну с шарами:

Урна — это все гены, которые мы в принципе могли обнаружить в нашем эксперименте (Background Universe, или фоновое множество). Пусть в урне $N = 20000$ шаров (генов).
Мы знаем, что путь «Гликолиз» состоит из $K = 100$ генов. Это красные шары в нашей урне. Остальные 19900 — белые.
В результате эксперимента алгоритм DESeq2 выдал нам список из $n = 500$ дифференциально экспрессирующихся генов. Мы достаем из урны 500 шаров.
Среди вытащенных шаров оказалось $k = 15$ красных (генов гликолиза).

Вероятность $P$ получить ровно $k$ успехов в выборке размера $n$ из генеральной совокупности $N$ , содержащей $K$ успешных элементов, описывается формулой гипергеометрического распределения:

P(X = k) = \frac{\binom{K}{k} \binom{N - K}{n - k}}{\binom{N}{n}}

Где $\binom{a}{b}$ — биномиальный коэффициент (число сочетаний). Нас интересует вероятность вытащить 15 или более красных шаров случайным образом. Если эта кумулятивная вероятность (p-value) крайне мала (например, $p < 0.01$ ), мы отвергаем нулевую гипотезу о случайном совпадении и заявляем, что путь «Гликолиз» статистически значимо обогащен (over-represented) в нашем списке DEGs.

Критическая ошибка ORA: неправильный выбор фона

Самая частая методологическая ошибка при проведении ORA, ломающая всю статистику — использование всего референсного генома в качестве фонового множества (Background Universe, параметр $N$ ).

Представьте, что вы анализируете транскриптом нейронов. В нейронах из-за эпигенетического сайленсинга в принципе не экспрессируются гены, специфичные для печени или мышц. Если вы используете как фон все 25000 генов человека, алгоритм будет думать: «Удивительно, в выборке совсем нет генов печени, зато так много генов нервной системы!». Вы получите мощное ложное обогащение базовых нейрональных путей, даже если ваш эксперимент изучал лишь тонкий механизм ответа на легкий стресс.

Правильный фон для ORA — это только те гены, которые имели шанс попасть в список DEGs в вашем конкретном эксперименте. Технически это гены, которые прошли первоначальный фильтр низких каунтов (например, гены, у которых baseMean > 0 в DESeq2). Если после фильтрации в матрице осталось 14000 экспрессирующихся генов, именно они, а не 25000, составляют $N$ в гипергеометрическом тесте.

Ограничения ORA

Несмотря на популярность и простоту реализации в пакетах вроде clusterProfiler или веб-инструментах типа DAVID, ORA имеет фундаментальные недостатки:

Зависимость от жестких порогов. Чтобы получить список DEGs, мы обязаны выбрать субъективные пороги. Ген с $FDR = 0.049$ попадет в анализ, а ген с $FDR = 0.051$ будет выброшен, хотя биологически между ними нет никакой разницы.
Игнорирование силы эффекта. ORA относится к списку DEGs как к плоскому бинарному множеству (ген либо есть в списке, либо нет). Ген, экспрессия которого изменилась в 100 раз, имеет ровно такой же вес при пересечении множеств, как ген, изменившийся в 1.5 раза.
Потеря слабых, но скоординированных сигналов. Если в метаболическом пути экспрессия всех 50 ферментов синхронно выросла на 15% (что недостаточно для преодоления порога FDR для каждого отдельного гена, но биологически приведет к мощному усилению синтеза метаболита), ORA вообще не увидит этот путь. Ни один ген не попадет в выборку, и путь будет пропущен.

Gene Set Enrichment Analysis (GSEA): анализ без порогов

Метод GSEA был разработан исследователями из Broad Institute именно для того, чтобы преодолеть ограничения ORA. Главная философия GSEA: не нужно отсекать гены искусственными порогами, нужно анализировать весь транскриптом целиком, оценивая сдвиги в распределении.

Шаг 1: Ранжирование генома

Вместо того чтобы делить гены на «значимые» и «незначимые», GSEA требует выстроить все обнаруженные в эксперименте гены (обычно 12000–18000) в единый непрерывный ранжированный список $L$ : от самых сильно активируемых до самых сильно подавляемых.

Какую метрику использовать для ранжирования? Использование только Log2 Fold Change (LFC) ошибочно: ген с огромным LFC, но единичными прочтениями (высоким шумом и высоким p-value) окажется на вершине списка, сломав анализ. Использование только p-value тоже не подходит, так как малые p-value бывают и у растущих, и у падающих генов (мы потеряем направление).

Золотой стандарт биоинформатики — комбинированная метрика $R$ , учитывающая и направление изменения, и его статистическую достоверность:

R = \text{sign}(LFC) \times -\log_{10}(p)

Где $\text{sign}(LFC)$ возвращает $1$ для генов с растущей экспрессией и $-1$ для генов с падающей. Функция $-\log_{10}(p)$ превращает маленькие p-value в большие положительные числа. Например, если ген имеет $LFC = -2.5$ и $p = 10^{-6}$ , его метрика $R = -1 \times 6 = -6$ . Он отправится в самый низ списка. В середине списка сгрудятся гены с $p \approx 1$ (их метрика $R \approx 0$ ) — это биологический шум, не отреагировавший на стимул.

Шаг 2: Алгоритм Running Sum

Имея отранжированный список всех генов $L$ , алгоритм GSEA берет конкретный биологический путь (набор генов $S$ , например, 80 генов из пути апоптоза) и начинает идти по нашему списку сверху вниз, вычисляя текущую сумму (Running Sum).

Правила шага:

Если алгоритм встречает ген, который принадлежит пути $S$ , сумма увеличивается (шаг вверх).
Если алгоритм встречает ген, который не принадлежит пути $S$ , сумма уменьшается (шаг вниз).

Величина шага вверх не одинакова — она взвешена по метрике ранжирования. Чем ближе ген к краям списка (где изменения сильнее), тем больший прирост он дает сумме. Шаг вниз рассчитывается так, чтобы при прохождении всего списка до конца сумма всегда возвращалась к нулю.

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

Шаг 3: Оценка результатов (ES, NES и Leading Edge)

Максимальное отклонение кривой от нуля называется Enrichment Score (ES). Однако сравнивать сырой ES между путями нельзя: путь из 300 генов математически имеет шанс набрать больший пик, чем путь из 20 генов, просто за счет количества шагов. Чтобы нивелировать размер пути, ES нормализуют путем перестановок (пермутаций), получая Normalized Enrichment Score (NES).

Положительный NES ( $> 1.5$ ) означает, что путь активирован (сдвинут в топ списка), отрицательный ( $< -1.5$ ) — подавлен (сдвинут в конец списка). Достоверность NES оценивается через FDR; обычно значимыми считаются пути с $FDR < 0.05$ или $0.25$ (в зависимости от строгости дизайна).

Важнейший практический выход GSEA — это Leading Edge subset (передовой край). Это подмножество генов исследуемого пути, которые встретились алгоритму до достижения максимального пика ES. Именно эти гены внесли основной вклад в сигнал обогащения. Если путь признан статистически значимым, биологу следует извлечь именно гены Leading Edge и построить по ним тепловую карту (heatmap) экспрессии — они являются главными молекулярными драйверами наблюдаемого фенотипа.

Сравнение подходов: что выбрать?

На практике ORA и GSEA отвечают на немного разные вопросы. В пайплайнах часто запускают оба метода: ORA для быстрой оценки самых сильных хитов, а GSEA — для поиска системных сдвигов.

Характеристика	Over-Representation Analysis (ORA)	Gene Set Enrichment Analysis (GSEA)
Входные данные	Короткий список значимых DEGs	Весь отранжированный транскриптом
Требование порогов	Да (жесткие отсечки по FDR и LFC)	Нет (используется непрерывная метрика)
Чувствительность	Низкая к слабым, но системным изменениям	Высокая к скоординированным слабым сдвигам
Устойчивость к шуму	Высокая (шумные гены отсекаются на этапе DESeq2)	Зависит от правильного выбора метрики ранжирования
Главный риск	Ошибка выбора фонового множества (Background)	Использование сырого LFC для ранжирования

Борьба с избыточностью результатов

Независимо от того, используете вы ORA или GSEA, при анализе по базе Gene Ontology вы неминуемо столкнетесь с проблемой избыточности. В результатах может оказаться 50 статистически значимых терминов, из которых 15 будут вариациями одного и того же процесса (например, "immune response", "innate immune response", "defense response to bacterium", "response to external biotic stimulus").

Это происходит из-за структуры DAG (правило истинного пути). Чтобы очистить результаты и оставить только смысловые кластеры, применяют алгоритмы редукции размерности на основе семантического сходства. Инструменты вроде пакета simplifyEnrichment (в R) или веб-сервера REVIGO вычисляют матрицу дистанций между терминами GO (на основе доли общих генов между ними) и схлопывают похожие термины в один репрезентативный кластер.

Это критически важный шаг перед визуализацией. Вместо нечитаемого barplot с сотней строк вы строите emapplot (Enrichment Map) или чистый dotplot из 5–10 главных биологических тем, затронутых в вашем эксперименте.

Основы single-cell RNA-Seq: технологические платформы, захват клеток и специфическая предобработка данных

Анализ транскриптома опухолевого биоптата классическим bulk RNA-Seq может показать умеренную, клинически незначимую экспрессию гена резистентности к химиотерапии. Врач назначает стандартный протокол, но через полгода происходит рецидив. Причина кроется в том, что опухоль — это не гомогенная масса. Если ген резистентности экстремально высоко экспрессируется лишь в 2% раковых стволовых клеток, а остальные 98% его не транскрибируют, усредненный сигнал bulk RNA-Seq размажет этот пик до уровня фонового шума. Получается «средняя температура по больнице», из-за которой упускается критически важная миноритарная субпопуляция. Технологии секвенирования РНК одиночных клеток (scRNA-Seq) решают эту проблему, превращая один биологический образец в тысячи независимых наблюдений, где каждая клетка становится отдельной точкой данных.

Технологические платформы: от лунки к капле

Чтобы прочитать транскриптом отдельной клетки, ее нужно физически изолировать от соседей, разрушить клеточную мембрану (лизис), захватить молекулы мРНК и провести обратную транскрипцию. Исторически и технологически сформировались два принципиально разных подхода к изоляции, каждый из которых определяет границы применимости метода.

Планшетные методы (Smart-seq2 / Smart-seq3)

В основе планшетных методов лежит физическая сортировка клеток. Клетки пропускаются через проточный цитофлуориметр (FACS) и по одной распределяются в лунки 96- или 384-луночного планшета. В каждой лунке независимо протекают реакции лизиса и подготовки библиотеки. В методах семейства Smart-seq используется обратная транскриптаза со свойством смены матрицы (template switching), что позволяет синтезировать полноразмерную кДНК.

Главное преимущество этого метода — полноразмерное секвенирование транскрипта (full-length). Поскольку мы не привязаны к одному концу молекулы, мы можем анализировать альтернативный сплайсинг, находить точечные мутации по всей длине гена (например, выявлять конкретные мутации в онкогене BRAF) и различать изоформы белков.

Недостатки обусловлены физическими ограничениями формата: низкая пропускная способность (сотни клеток на эксперимент) и высокая стоимость подготовки одной клетки (около 10–20 долл. США). Кроме того, работа с тысячами лунок требует сложной роботизированной станции для дозирования жидкостей.

Капельная микрофлюидика (10x Genomics, Drop-seq)

Современный стандарт высокопроизводительного scRNA-Seq опирается на микрофлюидные чипы. В чипе пересекаются три потока: суспензия одиночных клеток, суспензия микросфер (beads), покрытых олигонуклеотидами, и гидрофобное масло. На перекрестке потоков формируется эмульсия — микроскопические капли воды в масле (GEMs — Gel Bead-in-Emulsion). Идеальная капля содержит ровно одну клетку и ровно одну микросферу. Такая капля работает как изолированный нанореактор объемом в несколько пиколитра.

Капельные методы обладают колоссальной пропускной способностью (от 10 000 до 1 000 000 клеток в одном эксперименте) и низкой стоимостью (около 0.05–0.10 долл. за клетку). Однако из-за особенностей химии захвата они позволяют секвенировать только один конец транскрипта — обычно 3'-конец (реже 5'-конец). Молекула мРНК захватывается за поли-А хвост, и секвенируется лишь короткий фрагмент рядом с ним, что делает невозможным полноценный анализ сплайсинга.

Характеристика	Smart-seq2 (Планшеты)	10x Genomics (Микрофлюидика)
Пропускная способность	Сотни клеток	Десятки тысяч клеток
Покрытие транскрипта	Полноразмерное	Только 3' или 5' конец
Чувствительность (генов/клетку)	Высокая (3000–8000)	Средняя (1000–3000)
Анализ изоформ и мутаций	Да	Нет (только для концов)
Стоимость на 1 клетку	Высокая	Низкая

Математика инкапсуляции и проблема дуплетов

При загрузке клеток в микрофлюидный чип возникает статистическая проблема. Клетки распределяются по каплям случайно, и этот процесс математически описывается распределением Пуассона.

Вероятность $P(k)$ того, что в каплю попадет ровно $k$ клеток, вычисляется по формуле:

P(k) = \frac{\lambda^k e^{-\lambda}}{k!}

Где $\lambda$ — среднее ожидаемое количество клеток на одну каплю, а $e$ — основание натурального логарифма.

Если попытаться загрузить суспензию плотно, чтобы $\lambda = 1$ (в среднем одна клетка на каплю), мы получим:

Пустые капли ( $k=0$ ): $P(0) \approx 36.8\%$
Одиночные клетки ( $k=1$ ): $P(1) \approx 36.8\%$
Две и более клеток ( $k \geq 2$ ): $P(\geq 2) \approx 26.4\%$

Более четверти всех капель будут содержать несколько клеток — это дуплеты (doublets) или мультиплеты. Дуплеты катастрофичны для анализа: если в одну каплю попадут Т-лимфоцит и макрофаг, их РНК смешается. Биоинформатический алгоритм позже идентифицирует эту каплю как гибридный тип клетки, экспрессирующий маркеры и лимфоцитов, и макрофагов, которого не существует в природе.

Чтобы минимизировать количество дуплетов, клетки загружают в сильно разбавленном виде, искусственно снижая $\lambda$ . Например, при $\lambda = 0.1$ :

Пустые капли: $P(0) \approx 90.5\%$
Одиночные клетки: $P(1) \approx 9.0\%$
Дуплеты: $P(\geq 2) \approx 0.5\%$

Среди капель, содержащих хотя бы одну клетку, доля дуплетов составит около $5\%$ . Это осознанная плата за чистоту эксперимента: прибор генерирует миллионы пустых капель, чтобы те редкие капли, в которых есть клетки, содержали ровно одну.

Баркодирование: идентификация клеток и молекул

В капельной микрофлюидике все капли после лизиса и первичной реакции сливаются вместе (эмульсия разрушается), и дальнейшая подготовка библиотеки, а также само секвенирование происходят в одной общей пробирке. Чтобы после секвенирования понять, откуда взялся каждый короткий фрагмент ДНК (рид), используется двойная система молекулярных меток, расположенных на микросфере.

Cell Barcode (Клеточный баркод)

На одной микросфере закреплены миллионы зондов, и все они имеют одинаковую последовательность длиной 16 нуклеотидов — клеточный баркод (CB). Микросфер в коммерческом наборе миллионы, и каждая несет свой уникальный CB. Когда клетка лизируется внутри капли, все ее молекулы мРНК гибридизуются с зондами одной микросферы и получают одинаковый CB. Таким образом, CB отвечает на вопрос: «Из какой клетки пришел этот транскрипт?»

Unique Molecular Identifier (UMI)

Вторая метка — UMI (уникальный молекулярный идентификатор) — решает фундаментальную проблему количественной оценки. Перед секвенированием библиотека многократно амплифицируется с помощью полимеразной цепной реакции (ПЦР). Эффективность ПЦР зависит от длины и GC-состава фрагмента. Короткий ген с оптимальным GC-составом может дать 100 ПЦР-копий с одной исходной молекулы, а длинный ген со сложной вторичной структурой — всего 2 копии. Если просто посчитать итоговые риды, мы измерим эффективность ПЦР, а не реальную экспрессию гена в клетке.

UMI — это случайная последовательность из 10–12 нуклеотидов. В отличие от CB, который одинаков для всей микросферы, UMI уникален для каждого отдельного зонда на этой микросфере.

Когда оригинальная молекула мРНК связывается с зондом, она получает случайный UMI. При последующей ПЦР копируется весь конструкт (транскрипт + CB + UMI). После секвенирования биоинформатический пайплайн группирует риды не только по гену и клеточному баркоду, но и по UMI. Если пайплайн видит 100 ридов, картированных на ген GAPDH, с одинаковым клеточным баркодом и одинаковым UMI, он понимает, что это ПЦР-клоны одной и той же исходной молекулы. Эти 100 ридов «схлопываются» (collapsing) в 1 каунт. UMI отвечает на вопрос: «Сколько уникальных молекул мРНК этого гена было в клетке до ПЦР?»

Длина UMI в 12 нуклеотидов дает $4^{12} \approx 16.7 \times 10^6$ уникальных комбинаций. Риск того, что две разные молекулы одного и того же гена в одной клетке случайно свяжутся с зондами, имеющими одинаковый UMI (UMI collision), математически пренебрежим.

Специфическая предобработка данных

Геометрия чтения в капельном scRNA-Seq отличается от bulk-подходов. Секвенатор работает в асимметричном парноконцевом режиме. Read 1 (короткий, обычно 28 bp) содержит только техническую информацию: 16 bp клеточного баркода и 12 bp UMI. Read 2 (длинный, обычно 90 bp) содержит биологическую последовательность кДНК.

Программы первичной обработки (например, Cell Ranger от 10x Genomics или STARsolo) выполняют следующий пайплайн:

Выравнивание Read 2 на референсный геном (учитывая экзон-интронную структуру с помощью сплайс-ориентированных алгоритмов).
Присвоение каждому успешному выравниванию тегов CB и UMI из соответствующего Read 1.
Коррекция ошибок секвенирования в баркодах. Алгоритм сравнивает прочитанные CB с белым списком (whitelist) известных баркодов, заложенных производителем микросфер. Если допущена одна ошибка (расстояние Хэмминга равно 1), баркод исправляется.
Подсчет уникальных UMI для каждого гена внутри каждого клеточного баркода.

Результатом работы является матрица экспрессии (Count Matrix), где строки — это гены, а столбцы — клеточные баркоды. В отличие от bulk RNA-Seq, эта матрица экстремально разреженная (sparse): 80–90% ее значений составляют нули.

Разреженность возникает по двум причинам. Во-первых, клетка физически не экспрессирует все 20 000 генов генома одновременно. Во-вторых, из-за ограниченной эффективности захвата (capture efficiency около 10–20%) многие реально присутствующие транскрипты просто не попадают в библиотеку. Это явление называется эффектом drop-out. Из-за огромного количества нулей матрицы сохраняют в специальных форматах (например, Market Exchange Format, MTX), где записываются только ненулевые значения и их координаты, что экономит гигабайты оперативной памяти.

Идентификация реальных клеток (Cell Calling)

После генерации матрицы возникает проблема: столбцов (клеточных баркодов) в ней сотни тысяч, хотя в прибор загружалось всего 10 000 клеток. Откуда взялись остальные?

В суспензии всегда присутствует фоновая РНК (ambient RNA, или «суп»). Она высвобождается из клеток, которые погибли и разрушились еще на этапе диссоциации ткани до загрузки в чип. Эта свободная РНК плавает в растворе и неизбежно попадает в те самые «пустые» капли, в которых есть микросфера, но нет целой клетки. Фоновая РНК гибридизуется, получает баркоды и секвенируется.

Необходимо отделить баркоды, соответствующие реальным интактным клеткам, от баркодов пустых капель с фоновой РНК. Классический инструмент визуализации этого процесса — Knee plot (график колена).

Баркоды ранжируются по убыванию общего числа обнаруженных UMI и строятся на графике с двойной логарифмической шкалой.

На графике формируется характерный обрыв (колено). Баркоды слева от обрыва имеют тысячи и десятки тысяч UMI — это реальные клетки. Баркоды справа имеют 10–100 UMI — это пустые капли. Установка жесткого порога по количеству UMI (например, отсечение всех баркодов с $< 500$ UMI) работает для гомогенных клеточных линий.

Однако в реальных тканях присутствуют клетки разного размера. Активированный макрофаг огромен и содержит много РНК, а покоящийся T-лимфоцит мал и содержит мало РНК. Жесткий порог может ошибочно отфильтровать реальные мелкие клетки, приняв их за пустые капли.

Для решения этой проблемы применяются статистические алгоритмы, такие как EmptyDrops. Алгоритм работает в два этапа:

Берет пул баркодов с экстремально низким числом UMI (например, $< 100$ ), которые гарантированно являются пустыми каплями, и строит вероятностную модель (распределение Дирихле-мультиномиальное) профиля экспрессии фонового «супа».
Анализирует баркоды в спорной зоне (например, от 100 до 500 UMI). Если профиль экспрессии спорного баркода статистически неотличим от фонового супа, баркод удаляется. Если же профиль значимо отличается (например, содержит специфические маркеры лимфоцитов, которых мало в фоне), алгоритм сохраняет этот баркод как реальную клетку с низким содержанием РНК.

Контроль качества на уровне клеток (QC)

Идентификация баркодов — это лишь первый шаг. Даже если баркод содержит реальную клетку, она может быть непригодна для анализа. Перед нормализацией и кластеризацией матрицу подвергают жесткой фильтрации по трем основным метрикам.

Количество детектированных генов (nFeature_RNA). Клетки с аномально низким числом генов (например, $< 200$ ) обычно представляют собой сильно поврежденные клетки или артефакты секвенирования, где прочиталась лишь малая часть транскриптома. Клетки с аномально высоким числом генов часто являются скрытыми дуплетами, которые проскользнули через микрофлюидную систему (так как две слившиеся клетки разных типов дадут более широкое разнообразие генов, чем одна).
Общее число молекул / UMI (nCount_RNA). Метрика, тесно коррелирующая с числом генов. Слишком высокие значения также сигнализируют о потенциальных дуплетах.
Процент митохондриальных ридов (percent.mt). Это важнейший маркер клеточного стресса и апоптоза. Когда клетка повреждается в процессе диссоциации ткани, ее плазматическая мембрана рвется первой. Цитоплазматическая мРНК вытекает наружу (пополняя фоновый суп). Однако митохондрии имеют прочную двойную мембрану и остаются внутри клеточного каркаса. В результате в поврежденной клетке доля митохондриальных транскриптов искусственно взлетает с нормальных 1–5% до 20–50%. Такие «умирающие» клетки необходимо удалять, иначе их стрессовый транскриптомный профиль сформирует ложные кластеры при дальнейшем анализе.

Важный нюанс: порог для percent.mt не является универсальным. Для мононуклеарных клеток периферической крови (PBMC) нормой считается $< 5\%$ . Но для кардиомиоцитов (клеток сердечной мышцы), которые требуют огромного количества энергии и буквально набиты митохондриями, нормальным может быть показатель в 20–30%. Порог всегда подбирается индивидуально под тип ткани.

Процесс перехода от сырых FASTQ файлов к чистой, отфильтрованной от артефактов, пустых капель и мертвых клеток матрице экспрессии — фундамент scRNA-Seq исследования. Только убедившись, что каждая колонка в матрице представляет собой качественную одиночную клетку, мы можем переходить к математическому снижению размерности и поиску биологического смысла.

Математическое снижение размерности и алгоритмы кластеризации в анализе единичных клеток

В пространстве транскриптома млекопитающего около 20 000 измерений. Каждая ось — это уровень экспрессии отдельного гена, а каждая клетка в биологическом образце — точка, висящая в этом гиперобъеме. При попытке измерить классическое евклидово расстояние между любыми двумя случайно выбранными клетками в таком пространстве обнаруживается парадокс: все точки оказываются примерно на одинаковом расстоянии друг от друга. Это математическое явление называется «проклятием размерности». В высокоразмерных пространствах весь объем концентрируется в тонком слое у поверхности, алгоритмы машинного обучения теряют способность отличать близкие объекты от далеких, а вычислительная сложность матричных операций возрастает экспоненциально.

Чтобы найти биологический смысл — выделить уникальные типы клеток, отследить траектории их развития и найти редкие субпопуляции — необходимо сжать эти 20 000 измерений до нескольких десятков, а затем и до двух, не потеряв при этом истинные биологические связи. Этот процесс требует последовательного применения статистических фильтров, линейной алгебры и топологических алгоритмов.

Отбор признаков: фильтрация шума и поиск вариабельности

Первый шаг к снижению размерности делается до применения сложной тензорной математики и базируется на биологической логике. Из 20 000 генов значительная часть экспрессируется на стабильном базовом уровне во всех клетках. Это так называемые гены домашнего хозяйства (housekeeping genes), такие как ACTB (бета-актин) или GAPDH. Другая огромная часть генов в конкретной ткани не экспрессируется вообще (например, гены нейрональных рецепторов в клетках печени). Обе эти группы создают вычислительный шум и не помогают алгоритму отличить, например, Т-хелпер от цитотоксического Т-лимфоцита.

Для анализа отбирают высоко вариабельные гены (Highly Variable Genes, HVG). В данных RNA-Seq (особенно single-cell, где данные разрежены и содержат множество нулей, или drop-outs) дисперсия гена напрямую зависит от его среднего уровня экспрессии: чем выше средний каунт транскриптов, тем выше дисперсия. Это свойство отрицательного биномиального распределения, которым описываются данные секвенирования. Если отбирать гены просто по максимальной абсолютной дисперсии, в топ попадут только самые сильно экспрессирующиеся гены (например, рибосомальные белки), а критически важные, но редкие транскрипционные факторы (такие как FOXP3, определяющий регуляторные Т-клетки) будут проигнорированы.

Алгоритмы биоинформатических пакетов (функция FindVariableFeatures в Seurat или sc.pp.highly_variable_genes в Scanpy) моделируют ожидаемую зависимость между средним значением и дисперсией. Для этого применяется локальная полиномиальная регрессия (LOESS). Затем для каждого гена вычисляется стандартизированная дисперсия — метрика, показывающая, насколько реальная вариабельность гена превышает ожидаемую для его базового уровня экспрессии.

В результате исходная матрица размерностью $20000 \times N$ (где $N$ — число клеток) сокращается до матрицы $2000 \times N$ или $3000 \times N$ . Это первое, биологически обоснованное снижение размерности, которое убирает фоновый шум и оставляет только те векторы, которые содержат информацию о клеточной гетерогенности.

Метод главных компонент (PCA) как фундамент анализа

Даже 2000 измерений — это слишком много для точного вычисления расстояний и построения графов. Следующим шагом применяется метод главных компонент (Principal Component Analysis, PCA). Это строгий линейный алгоритм, который ищет в данных новые ортогональные оси (главные компоненты), вдоль которых дисперсия данных максимальна.

Математически PCA сводится к вычислению ковариационной матрицы признаков и нахождению ее собственных векторов и собственных значений. Если $X$ — центрированная матрица экспрессии (где по строкам расположены клетки, а по столбцам — отобранные вариабельные гены), то ковариационная матрица $C$ вычисляется как:

C = \frac{1}{n-1} X^T X

где:

$C$ — ковариационная матрица генов, отражающая их совместную изменчивость.
$n$ — количество клеток в эксперименте.
$X$ — матрица экспрессии, где среднее значение каждого гена сдвинуто к нулю.
$X^T$ — транспонированная матрица экспрессии.

Собственные векторы матрицы $C$ задают направления новых осей (PC1, PC2, PC3 и так далее), а собственные значения показывают, какую долю общей дисперсии объясняет каждая конкретная компонента. Первая компонента (PC1) всегда строится так, чтобы объяснить максимальный объем вариации в данных. На практике в scRNA-Seq PC1 часто разделяет клетки по глобальным признакам, например, отделяет иммунные клетки от эпителиальных. PC2 объясняет максимальный из оставшегося объема дисперсии, при строгом условии ортогональности (перпендикулярности) к PC1.

В анализе единичных клеток PCA играет критическую роль мощного фильтра шумоподавления. Истинный биологический сигнал — различия между типами и подтипами клеток — обычно улавливается первыми 15–50 компонентами. Все последующие компоненты (от 51 до 2000) содержат преимущественно технический шум эксперимента: случайные флуктуации эффективности обратной транскрипции, вариации амплификации и ошибки секвенирования. Оставляя только первые 30 компонент, мы проецируем клетки в 30-мерное пространство. Именно в этом плотном, очищенном от шума пространстве в дальнейшем будут вычисляться расстояния между клетками для графовой кластеризации.

Выбор количества главных компонент — задача, требующая баланса. Если взять слишком мало компонент (например, 5), алгоритм «недообучится» на биологии: редкие клеточные популяции, такие как дендритные клетки в крови, чьи уникальные маркеры формируют, скажем, 12-ю и 15-ю компоненты, сольются с макрофагами. Если взять слишком много (например, 100), алгоритм начнет искать кластеры в техническом шуме. На практике используют график «каменистой осыпи» (Elbow plot), где по оси X отложен номер компоненты, а по оси Y — доля объясненной дисперсии. Точка перегиба («локоть»), после которой график становится пологим и стремится к нулю, указывает на оптимальное число компонент.

Нелинейное снижение размерности: топология и визуализация

Биологические процессы редко развиваются по прямым линиям. Дифференцировка стволовой клетки крови в зрелый эритроцит — это непрерывный каскад активации и репрессии тысяч генов. В высокоразмерном пространстве такие процессы формируют сложные, изогнутые структуры — нелинейные многообразия (manifolds).

PCA, будучи линейным методом, проецирует данные как тени на плоскую стену. Если многообразие имеет форму скрученного рулета (классический датасет Swiss Roll), линейная проекция наложит разные слои этого рулета друг на друга. Клетки, находящиеся на разных стадиях развития, на графике PC1 vs PC2 могут оказаться в одной точке.

Чтобы аккуратно развернуть такие структуры на 2D-плоскость для визуализации человеком, применяются алгоритмы manifold learning. Самые востребованные в биоинформатике — t-SNE и UMAP.

t-SNE и решение проблемы скученности

Алгоритм t-SNE (t-distributed Stochastic Neighbor Embedding) совершил революцию в визуализации транскриптомов. Его фундаментальная идея — полный отказ от сохранения абсолютных дистанций между точками. Вместо этого алгоритм фокусируется на вероятности того, что две клетки являются ближайшими соседями.

В исходном 30-мерном PCA-пространстве t-SNE измеряет сходство между клетками с помощью нормального (гауссовского) распределения. Однако при попытке перенести эти вероятности в 2D-пространство возникает «проблема скученности» (crowding problem). В 30-мерном пространстве объем сферической окрестности вокруг точки растет колоссальными темпами. На плоскости места катастрофически не хватает. Если попытаться перенести все многомерные расстояния пропорционально, точки на плоскости сольются в единое плотное пятно в центре координат.

Чтобы компенсировать нехватку пространства, t-SNE использует для 2D-проекции распределение Стьюдента с одной степенью свободы (распределение Коши). В отличие от Гауссианы, оно имеет «тяжелые хвосты». Математически это означает, что для сохранения той же вероятности соседства, клетки, находящиеся на средних дистанциях в многомерном пространстве, на 2D-плоскости должны быть оттолкнуты друг от друга значительно дальше. Благодаря этому механизму отталкивания t-SNE блестяще разрывает разные типы клеток на четкие, визуально изолированные островки.

Ключевой гиперпараметр t-SNE — перплексия (perplexity). Практически он означает ожидаемое количество ближайших соседей у каждой клетки, на которые алгоритм обращает внимание.

Низкая перплексия (например, 5) заставляет алгоритм фокусироваться только на самых близких соседях. Локальная структура доминирует, что часто приводит к распаду единой биологической популяции на множество мелких бессмысленных кластеров. Слишком высокая перплексия (например, 500) заставляет алгоритм учитывать слишком много дальних связей, из-за чего разные типы клеток слипаются. Оптимальное значение обычно лежит в диапазоне от 30 до 50.

Главный недостаток t-SNE — полное разрушение глобальной топологии. Расстояние между двумя разными кластерами на графике t-SNE не значит ничего. Если кластер B-клеток визуально находится рядом с кластером фибробластов, это не делает их биологически родственными — это артефакт математической развертки.

UMAP: алгебраическая топология и глобальная структура

UMAP (Uniform Manifold Approximation and Projection) опирается на риманову геометрию. В отличие от t-SNE, который оперирует исключительно вероятностями, UMAP строит симплициальные комплексы (графы, включающие узлы, ребра и заполненные объемы), аппроксимируя геометрическую форму данных в многомерном пространстве, а затем пытается собрать максимально похожий комплекс на плоскости.

UMAP вытеснил t-SNE и стал стандартом де-факто по двум причинам. Во-первых, вычислительная оптимизация позволяет ему работать на порядки быстрее, что критично для современных датасетов на сотни тысяч и миллионы клеток. Во-вторых, за счет инициализации графа через спектральное вложение (Spectral embedding), UMAP гораздо лучше сохраняет глобальную структуру данных. Расстояния между кластерами обретают смысл: родственные субпопуляции (например, CD4+ и CD8+ Т-клетки) с высокой вероятностью окажутся рядом, образуя единый суперкластер, а далекие по происхождению клетки (нейроны и макрофаги) — на разных концах графика. Непрерывные процессы, такие как созревание клеток, образуют на UMAP красивые связные траектории.

Важное правило транскриптомного анализа: t-SNE и UMAP используются исключительно для визуализации. Никогда не запускайте алгоритмы кластеризации или анализ дифференциальной экспрессии генов, используя 2D-координаты UMAP. При экстремальном сжатии из 30 измерений в 2 неизбежно происходят искажения, наложения и потеря информации. Вся строгая математика должна выполняться в многомерном пространстве главных компонент (PCA).

Графовая кластеризация: от матриц к сетям

Имея чистое 30-мерное PCA-пространство, необходимо сгруппировать клетки в биологические типы. Классические алгоритмы машинного обучения, такие как K-means, здесь терпят неудачу. K-means математически предполагает, что кластеры имеют сферическую форму и примерно одинаковый размер. Клеточные же популяции часто образуют вытянутые траектории (клетки в процессе деления) или имеют колоссальную разницу в плотности (10 000 эритроцитов и 50 стволовых клеток). Кроме того, K-means требует заранее указать количество кластеров ( $K$ ), которое в исследовательском (exploratory) анализе новой ткани неизвестно.

Поэтому современным стандартом стала графовая кластеризация (Graph-based clustering), не требующая знания о количестве кластеров заранее.

Шаг 1: Построение KNN-графа

В пространстве PCA для каждой клетки вычисляются евклидовы расстояния до всех остальных клеток и находятся $k$ ее ближайших соседей (обычно $k = 20 \dots 30$ ). Каждая клетка становится узлом графа, а ненаправленные связи проводятся к ее соседям. На этом этапе граф получается зашумленным: случайная клетка с техническими артефактами может оказаться чьим-то соседом просто из-за отсутствия других вариантов в пустом участке гиперобъема.

Шаг 2: Переход к SNN-графу (Shared Nearest Neighbor)

Чтобы сделать граф устойчивым к шуму, связи переоцениваются на основе общих соседей. Логика проста: если две клетки соединены в KNN-графе, но у них нет больше ни одного общего соседа, их связь случайна. Если же у них 15 общих соседей из 20 возможных, они явно принадлежат к одной плотной, биологически однородной популяции.

Вес связи между клетками вычисляется с помощью индекса Жаккара:

$J(A, B) = \frac{|A \cap B|}{|A \cup B|}$

где:

$J(A, B)$ — индекс Жаккара, мера сходства двух множеств (от 0 до 1).
$A$ — множество ближайших соседей первой клетки.
$B$ — множество ближайших соседей второй клетки.
$|A \cap B|$ — количество общих соседей (мощность пересечения множеств).
$|A \cup B|$ — общее количество уникальных соседей обеих клеток (мощность объединения множеств).

Например, если у клетки А 20 соседей, у клетки B 20 соседей, и 15 из них общие, то объединение составит $20 + 20 - 15 = 25$ уникальных клеток. Индекс Жаккара будет равен $15 / 25 = 0.6$ . Ребрам с низким индексом Жаккара присваивается вес 0 (связь обрывается). В результате получается взвешенный SNN-граф, где веса ребер отражают истинную топологическую близость транскриптомных профилей.

Шаг 3: Оптимизация модулярности (Louvain и Leiden)

Задачу кластеризации теперь можно сформулировать как классический поиск сообществ (communities) в графе. Алгоритм должен разрезать граф так, чтобы внутри групп связи были максимально плотными и тяжелыми, а между группами — редкими и легкими. Мерой качества такого разбиения служит модулярность ( $Q$ ).

Q = \frac{1}{2m} \sum_{i,j} \left[ A_{ij} - \frac{k_i k_j}{2m} \right] \delta(c_i, c_j)

где:

$Q$ — модулярность графа (значение от -1 до 1).
$m$ — сумма весов всех ребер в графе.
$A_{ij}$ — фактический вес ребра между узлами $i$ и $j$ .
$k_i$ и $k_j$ — суммы весов всех ребер, присоединенных к узлам $i$ и $j$ соответственно.
$\frac{k_i k_j}{2m}$ — ожидаемый вес ребра между $i$ и $j$ , если бы связи в графе распределялись абсолютно случайно.
$c_i, c_j$ — кластеры, к которым отнесены узлы $i$ и $j$ .
$\delta(c_i, c_j)$ — функция Кронекера: равна 1, если узлы в одном кластере ( $c_i = c_j$ ), и 0 в противном случае.

Смысл формулы: алгоритм суммирует разницу между реальным весом связи и тем весом, который ожидался бы при случайном перемешивании графа. Суммирование идет только для пар узлов внутри одного кластера (благодаря функции Кронекера, обнуляющей межкластерные связи). Чем выше $Q$ , тем сильнее связи внутри кластеров превышают случайные ожидания.

Алгоритм Louvain работает итеративно: сначала каждая клетка считается микро-кластером. Затем алгоритм пытается объединить соседние узлы, если это действие увеличивает глобальную модулярность $Q$ . Процесс повторяется, узлы сливаются в мета-узлы, пока $Q$ не достигнет локального максимума. У алгоритма Louvain есть известный математический изъян: при слиянии крупных мета-узлов он может создавать внутренне несвязные кластеры (клетки формально в одном кластере, но пути по графу между ними нет). Алгоритм Leiden решает эту проблему. На каждом шаге укрупнения графа Leiden позволяет временно разбивать уже сформированные сообщества, гарантируя, что итоговые кластеры будут строго связными. В современных пайплайнах (Seurat v4+, Scanpy) Leiden является методом по умолчанию.

Управление детализацией: параметр Resolution

Главный инструмент биолога при управлении графовой кластеризацией — параметр разрешения (resolution). Математически он внедряется в формулу модулярности как множитель перед штрафом за случайные связи (перед дробью $\frac{k_i k_j}{2m}$ ).

При низком разрешении (например, 0.2) алгоритм сильно штрафует за создание новых кластеров. В результате клетки объединяются в крупные макро-популяции: все Т-клетки сливаются в один кластер, все В-клетки — в другой. При высоком разрешении (1.0 и выше) штраф снижается, и алгоритм начинает дробить граф на мельчайшие субпопуляции. Т-клетки могут разделиться на наивные (CD4+ Naive), клетки памяти (Memory T), эффекторные и истощенные (Exhausted T-cells). Выбор оптимального разрешения не имеет правильного математического ответа — он всегда зависит от биологической задачи исследования.

Переход от сырых чтений секвенатора к осмысленным клеточным кластерам — это путь последовательного снижения размерности. Сначала отсекаются неинформативные гены (HVG), затем сжимаются линейные зависимости для удаления технического шума (PCA), и, наконец, используется топология данных для поиска плотных сообществ (SNN-граф и Leiden) с их последующей визуализацией (UMAP). Полученные кластеры — это строгие математические абстракции. На следующем этапе анализа предстоит вдохнуть в них биологический смысл, определив с помощью дифференциальной экспрессии, какие именно типы клеток скрываются за номерами кластеров.

Материалы для углубленного изучения математического аппарата

Для самостоятельного погружения в математические основы транскриптомики и алгоритмы машинного обучения, используемые в scRNA-Seq, рекомендуется обратиться к следующим ресурсам:

Фундаментальные учебники и книги:

Modern Statistics for Modern Biology (Susan Holmes, Wolfgang Huber) — классический англоязычный учебник, детально разбирающий применение PCA, кластеризации и статистических тестов на биологических данных. Доступен бесплатно онлайн.
Машинное обучение и анализ данных — материалы Школы анализа данных (ШАД) Яндекса. Отличный русскоязычный ресурс для глубокого понимания линейной алгебры, лежащей в основе PCA, и метрических алгоритмов кластеризации.

Руководства и документация (GitHub / Web): 3. Orchestrating Single-Cell Analysis with Bioconductor (OSCA) — исчерпывающая онлайн-книга (EN) от разработчиков пакетов на R. Содержит глубокий разбор математики отбора признаков (HVG) и графовой кластеризации. 4. Официальные туториалы пакета Scanpy (Python) на GitHub и ReadTheDocs. В разделе API подробно описана математика под капотом функций sc.tl.pca, sc.tl.umap и sc.tl.leiden. 5. Официальные виньетки пакета Seurat (R) от Satija Lab. Включают детальные объяснения выбора параметров (таких как resolution и perplexity) на реальных датасетах.

Оригинальные научные публикации (must-read для понимания алгоритмов): 6. UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction (Leland McInnes et al., 2018) — оригинальная статья, описывающая топологические основы UMAP. 7. From Louvain to Leiden: guaranteeing well-connected communities (Traag et al., 2019) — статья, математически доказывающая превосходство алгоритма Leiden над Louvain в графовой кластеризации.

Идентификация клеточных типов, поиск маркерных генов и автоматизированная аннотация популяций

После успешного снижения размерности и графовой кластеризации исследователь неизбежно сталкивается с экраном, на котором россыпь точек разбита на аккуратные, но совершенно безликие группы: «Кластер 0», «Кластер 1», «Кластер 2». Математика выполнила свою задачу, сгруппировав транскриптомно схожие клетки. Однако для биологии эти номера не значат ничего. В классическом транскриптомном анализе (bulk RNA-Seq) мы заранее знаем, что секвенируем: ткань печени сравнивается с тканью легкого, или опухоль сравнивается со здоровым контролем. В single-cell RNA-Seq (scRNA-Seq) парадигма переворачивается: мы сначала секвенируем «суп» из клеток, кластеризуем их вслепую, и лишь затем выясняем, с чем именно имеем дело.

Процесс перевода математических кластеров в биологические термины (Т-лимфоциты, астроциты, кардиомиоциты) называется аннотацией клеточных типов, и он требует совершенно иных статистических подходов, нежели анализ смешанных тканей.

Статистика поиска маркеров: почему bulk-методы здесь ломаются

В классическом bulk RNA-Seq стандартом де-факто для поиска дифференциально экспрессирующихся генов являются алгоритмы, основанные на отрицательном биномиальном распределении (например, DESeq2 или edgeR). Они моделируют ожидаемую дисперсию гена в зависимости от среднего уровня его экспрессии. Логично предположить, что для поиска маркеров кластера в single-cell данных можно применить тот же пайплайн: взять «Кластер 0» как опытную группу, а все остальные клетки объединить в контрольную (стратегия One-vs-All).

На практике эта стратегия сталкивается с тремя фундаментальными препятствиями scRNA-Seq:

Экстремальная разреженность данных (sparsity) и drop-out эффекты. В матрице scRNA-Seq от 70% до 90% значений — это нули. Отрицательное биномиальное распределение плохо справляется с таким избытком нулей, что приводит к неадекватной оценке дисперсии.
Колоссальный размер выборки. В bulk-эксперименте мы сравниваем, например, 3 биологических повторности с 3 контролями. В single-cell мы сравниваем 5000 клеток одного кластера с 15000 клеток других кластеров. При выборках такого объема ( $N \to \infty$ ) классические параметрические тесты становятся гиперчувствительными: малейшее, биологически абсолютно незначимое изменение экспрессии гена на доли процента получает астрономически малый $p$ -value (например, $10^{-250}$ ).
Вычислительная сложность. DESeq2 требует подгонки параметров модели для каждого гена. Расчет такой матрицы на десятки тысяч клеток может занять часы или даже дни.

Поэтому золотым стандартом для поиска маркерных генов в scRNA-Seq стал непараметрический критерий суммы рангов Уилкоксона (Wilcoxon rank-sum test), также известный как U-критерий Манна-Уитни. Его главное преимущество — работа не с абсолютными значениями каунтов, которые сильно искажены техническим шумом, а с их рангами.

Математическая логика критерия опирается на вычисление статистики $U$ . Для оценки того, экспрессируется ли ген выше в Кластере А по сравнению с Кластером Б (или всеми остальными клетками), алгоритм объединяет значения экспрессии этого гена из всех клеток, сортирует их по возрастанию и присваивает каждому значению ранг (от 1 до общего числа клеток). Затем ранги клеток, принадлежащих только Кластеру А, суммируются.

Статистика вычисляется по формуле:

U = R_A - \frac{n_A(n_A + 1)}{2}

Где $R_A$ — сумма рангов значений экспрессии в Кластере А, а $n_A$ — количество клеток в Кластере А. Элемент $\frac{n_A(n_A + 1)}{2}$ представляет собой минимально возможную сумму рангов для выборки размера $n_A$ (ситуация, когда все клетки Кластера А имеют самые низкие значения экспрессии).

Переход к рангам делает тест невероятно устойчивым к выбросам. Если из-за локальной ошибки ПЦР-амплификации одна клетка получила 5000 каунтов по гену, тогда как остальные имеют от 0 до 10, параметрический тест сместит среднее значение всей группы. В ранговом тесте эта клетка просто получит наивысший ранг (например, ранг 20000), не искажая общую статистику. Если ген является сильным маркером Кластера А, большинство его клеток займут верхние строчки в отсортированном списке, сумма рангов $R_A$ будет максимальной, и значение $U$ укажет на высокую значимость.

Анатомия идеального маркерного гена

Статистическая значимость ( $p$ -value) — лишь первый, самый грубый фильтр. В single-cell анализе ген может иметь $p$ -value $< 10^{-50}$ , но при этом быть абсолютно бесполезным для биологической аннотации. Идеальный маркерный ген должен обладать высокой дискриминационной способностью, которая оценивается через сочетание трех ключевых метрик.

Первая метрика — Log2 Fold Change (LFC). Она показывает силу биологического эффекта (во сколько раз экспрессия гена в кластере выше, чем вне его). Однако в scRNA-Seq LFC часто занижен из-за огромного количества нулей в фоновой группе.

Гораздо важнее две другие метрики, специфичные для анализа единичных клеток: $pct.1$ и $pct.2$ .

$pct.1$ — доля клеток внутри исследуемого кластера, в которых экспрессия данного гена строго больше нуля.
$pct.2$ — доля клеток во всех остальных кластерах (фоне), где этот ген детектирован.

Идеальный маркер имеет $pct.1 \approx 1.0$ (экспрессируется в 100% клеток кластера) и $pct.2 \approx 0.0$ (не экспрессируется нигде больше). На практике такие бинарные маркеры встречаются редко, и биологу приходится балансировать между чувствительностью и специфичностью.

Сравним три потенциальных маркера для гипотетического кластера Т-клеток:

Ген CD3D: $pct.1 = 0.95$ , $pct.2 = 0.40$ , $LFC = 1.5$ .
Ген FOXP3: $pct.1 = 0.15$ , $pct.2 = 0.01$ , $LFC = 3.0$ .
Ген ACTB: $pct.1 = 0.99$ , $pct.2 = 0.98$ , $LFC = 0.2$ .

Ген CD3D является отличным базовым маркером (пан-Т-клеточным): он охватывает почти весь кластер, но из-за биологического родства «размазан» и по другим популяциям (например, NK-клеткам). Ген FOXP3, напротив, крайне специфичен (почти не встречается вне кластера), но детектируется лишь в 15% клеток самого кластера. Это типичная картина для транскрипционных факторов: они биологически определяют тип клетки (регуляторные Т-клетки), но из-за низкой базовой экспрессии часто теряются в drop-out. Ген ACTB (бета-актин) экспрессируется везде. Несмотря на то, что он может быть статистически значимо выше в одном кластере, для аннотации он бесполезен.

Для формализации дискриминационной способности часто используют метрику AUC (Area Under the ROC Curve). Алгоритм пытается использовать экспрессию одного гена как бинарный классификатор: если экспрессия выше определенного порога — это клетка нашего кластера, если ниже — чужого.

$AUC = 0.5$ означает, что ген предсказывает принадлежность к кластеру не лучше подбрасывания монетки (как ген ACTB).
$AUC = 1.0$ означает идеальное разделение без ложноположительных и ложноотрицательных результатов. Гены с $AUC > 0.8$ обычно считаются надежными и устойчивыми маркерами.

Ручная аннотация и паттерны визуализации

Получив таблицу с топовыми маркерами для каждого кластера (отсортированную по LFC или AUC), исследователь приступает к сопоставлению найденных генов с известными биологическими базами данных и литературой.

Для визуальной оценки качества маркеров в статьях и отчетах используется специфический тип графика — Dot plot (пузырьковая диаграмма). Он решает проблему одновременного отображения двух параметров (уровня экспрессии и процента экспрессирующих клеток) для множества генов и множества кластеров.

В классическом Dot plot по оси X располагаются интересующие нас маркерные гены, по оси Y — номера кластеров. На пересечении рисуется круг. Размер круга строго привязан к параметру $pct.1$ (доле клеток кластера, экспрессирующих ген). Цвет круга (обычно градиент от светло-серого к ярко-красному или темно-синему) отражает средний уровень нормализованной экспрессии этого гена только среди тех клеток кластера, где он детектирован. Хорошо аннотированный Dot plot выглядит как четкая диагональ из крупных ярких кругов, где каждый клеточный тип имеет свой уникальный набор «горящих» маркеров.

Классический пример ручной аннотации — анализ мононуклеарных клеток периферической крови (PBMC). Если таблица маркеров для Кластера 2 выдает в топе гены CD79A, MS4A1 (кодирует белок CD20) и CD19, биолог переименовывает «Кластер 2» в «В-лимфоциты». Если Кластер 5 характеризуется высокой экспрессией LYZ, CD14 и S100A9, он аннотируется как «Моноциты».

Этот подход работает для любых тканей. При анализе коры головного мозга кластер с экспрессией SYT1 и SNAP25 будет определен как нейроны, а кластер с GFAP и AQP4 — как астроциты.

Сложности начинаются при попытке глубокой детализации. Т-клетки (CD3D+) могут разбиваться на десяток субкластеров. Различить наивные CD4+ Т-клетки и CD4+ Т-клетки памяти ручным перебором генов становится крайне сложно. Их транскриптомные профили различаются минимально, а классические поверхностные белки-маркеры, используемые в проточной цитометрии (CD45RA и CD45RO), являются изоформами сплайсинга одного и того же гена PTPRC. Стандартный 3'-scRNA-Seq протокол секвенирует только концы транскриптов и физически не способен различить эти изоформы. В таких случаях ручная аннотация упирается в технологический предел.

Автоматизированная аннотация: перенос ярлыков с референса

Когда ручной анализ становится неэффективным (например, при аннотации атласа на миллион клеток или поиске редких субпопуляций), применяются методы автоматизированной аннотации. Их фундаментальная идея — Reference-based mapping (картирование на эталон).

Суть метода: берется внешний, уже размеченный экспертами датасет (референс). Это может быть массивный bulk RNA-Seq очищенных клеточных популяций или качественный scRNA-Seq атлас здоровой ткани. Затем алгоритм математически сравнивает транскриптом каждой неизвестной клетки из нашего (query) датасета с профилями из референса.

Один из самых популярных алгоритмов этого класса — SingleR. Он работает на уровне единичных клеток, а не усредненных кластеров, что позволяет находить редкие подтипы клеток даже внутри смешанного графового кластера.

Механика SingleR построена на итеративном вычислении коэффициента ранговой корреляции Спирмена. Выбор ранговой корреляции не случаен: она устойчива к эффектам батча (batch effects) и разнице в глубине секвенирования между нашим датасетом и референсом.

На первом шаге алгоритм отбирает маркерные гены, которые вариабельны в самом референсном датасете.
Вычисляется корреляция Спирмена между профилем экспрессии неизвестной клетки и медианным профилем каждого клеточного типа в референсе.
Клетке присваивается предварительная метка того типа, с которым корреляция максимальна.
Затем SingleR выполняет fine-tuning (тонкую настройку). Если у клетки одинаково высокая корреляция с «Т-клетками памяти» и «Наивными Т-клетками», алгоритм отбрасывает все гены, кроме тех, которые дифференциально экспрессируются именно между этими двумя близкими типами, и пересчитывает корреляцию заново. Итерации продолжаются, пока не останется один однозначный победитель.

Альтернативный, более геометрический подход реализован в экосистеме Seurat (функции FindTransferAnchors и TransferData). Этот метод использует концепцию взаимных ближайших соседей (MNN — Mutual Nearest Neighbors). Алгоритм проецирует оба датасета в общее пространство сниженной размерности (обычно PCA) и ищет «якоря». Якорь — это пара клеток (одна из референса, другая из запроса), которые являются ближайшими соседями друг для друга в этом многомерном пространстве. Если клетка X из нашего опыта считает клетку Y из референса самой похожей на себя, а клетка Y отвечает взаимностью, между ними устанавливается связь. На основе сети таких якорей метки переносятся с референса на наш датасет с вычислением вероятности (prediction score).

Автоматическая аннотация невероятно ускоряет работу, но имеет критическую системную уязвимость: алгоритм не способен предсказать клеточный тип, которого изначально нет в референсе. Если в образце глиобластомы присутствует уникальная популяция опухоль-ассоциированных макрофагов, а референс построен на здоровом мозге, автоматика принудительно «натянет» на них наиболее похожий ярлык (например, обычная микроглия), скрыв от исследователя потенциальное открытие.

Анатомия «мусорных» и переходных кластеров

В реальном проекте после аннотации всегда остаются кластеры, маркерные профили которых выглядят как биологический абсурд. Умение интерпретировать такие аномалии и очищать от них данные — важнейший навык биоинформатика.

Первый тип аномалии — двойные маркеры (дуплеты). Если Кластер 7 одновременно экспрессирует мощные маркеры эпителиальных клеток (EPCAM, KRT18) и макрофагов (CD68, C1QA), это почти наверняка кластер нераспознанных дуплетов. Несмотря на работу алгоритмов фильтрации (таких как DoubletFinder или Scrublet) на ранних этапах, гетеротипические дуплеты — физическое слияние двух разных клеток в одной капле микрофлюидного чипа — часто выживают. На графике UMAP они формируют отдельный «мост» между двумя чистыми популяциями. Такие кластеры подлежат полному удалению из дальнейшего анализа.

Второй тип — стрессовые кластеры. Они не имеют четких маркеров клеточной идентичности (ни нейроны, ни фибробласты, ни Т-клетки), но в их профиле абсолютно доминируют митохондриальные гены, гены теплового шока (семейство HSP) и гены апоптоза. Это клетки, которые были критически повреждены ферментами в процессе диссоциации ткани, но не успели разрушиться до попадания в секвенатор. Их транскриптом отражает предсмертную агонию, а не реальную биологическую функцию.

Третий, самый сложный тип — переходные состояния. Иногда кластер не имеет собственных сильных маркеров, но обладает градиентом маркеров соседних кластеров. Например, в данных развивающегося костного мозга можно найти популяцию, плавно теряющую транскрипционные факторы стволовых клеток и одновременно набирающую маркеры эритроцитов. Это не техническая ошибка, а запечатленный в моменте процесс клеточной дифференцировки. Статичная кластеризация и жесткая аннотация (одна клетка = один тип) здесь пасуют. Для анализа таких непрерывных процессов требуются алгоритмы вывода траекторий (trajectory inference), которые будут рассмотрены на следующих этапах курса.

Аннотация клеточных типов — это момент, когда сухая статистика передает эстафету биологии. Ни один алгоритм не способен выдать абсолютную истину. Любая метка, присвоенная кластеру, будь то результат вдумчивого анализа Dot plot или работы SingleR — это лишь статистически обоснованная гипотеза, которая должна согласовываться с дизайном эксперимента и здравым смыслом исследователя.

Траекторный анализ и моделирование динамики генной экспрессии в процессах клеточной дифференцировки

Фундаментальный парадокс транскриптомики заключается в том, что для измерения уровня экспрессии генов в клетке мы обязаны эту клетку разрушить. Секвенирование РНК — это всегда посмертный снимок. Мы не можем взять одну конкретную стволовую клетку и наблюдать, как меняется ее транскриптом на протяжении пяти дней дифференцировки в нейрон. Однако биологические процессы, такие как эмбриогенез, иммунный ответ или регенерация тканей, разворачиваются во времени непрерывно. Если в наших руках есть только статические снимки десятков тысяч мертвых клеток, как мы можем реконструировать динамику их жизни?

В классическом bulk RNA-Seq эта задача неразрешима в принципе. Усредняя сигнал по миллионам асинхронно развивающихся клеток, мы получаем «транскриптомный шум», в котором смешаны ранние, промежуточные и поздние стадии. Переход к single-cell RNA-Seq (scRNA-Seq) позволил изолировать сигнал каждой клетки, однако проблема статического снимка осталась. Ответ на вызов кроется в математическом моделировании непрерывных многообразий и элегантном использовании биологических особенностей созревания матричной РНК.

Эргодическая гипотеза и концепция псевдовремени

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

Метафора марафона отлично иллюстрирует этот принцип. Моментальная фотография трассы фиксирует тысячи бегунов разного уровня подготовки. Кто-то только стартовал, кто-то преодолевает половину дистанции, а кто-то уже финиширует. Выстроив бегунов по их координатам на трассе, можно реконструировать весь процесс бега от старта до финиша, хотя фотография была сделана за долю секунды. Биологическим аналогом является сперматогенез: в семенниках одновременно присутствуют сперматогонии, сперматоциты разных порядков, сперматиды и зрелые сперматозоиды. Секвенировав эту ткань, мы получаем срез всех возможных состояний процесса.

В scRNA-Seq «трассой» выступает многомерное пространство генной экспрессии. Клетки, находящиеся на соседних этапах дифференцировки, имеют очень похожие транскриптомные профили. Выстраивая их в непрерывную цепочку на основе транскриптомного сходства, алгоритмы вычисляют одномерную координату — псевдовремя (pseudotime).

Псевдовремя — это не реальное физическое время в часах или сутках. Это мера транскриптомного прогресса клетки вдоль определенного биологического пути. Расстояние в псевдовремени между двумя клетками отражает объем транскрипционных изменений, которые необходимо совершить, чтобы перейти из одного состояния в другое.

Алгоритмы топологического вывода: как строятся траектории

Алгоритмы графовой кластеризации (Louvain, Leiden), реализованные в пакетах Seurat или Scanpy, разбивают клетки на дискретные группы. Это корректно для терминально дифференцированных клеток (например, зрелые B-клетки и Т-клетки в периферической крови). Но в развивающихся тканях границы между кластерами размыты: клетки образуют непрерывный континуум (manifold). Принудительное разбиение такого континуума на жесткие кластеры приводит к потере информации о переходных состояниях.

Для решения этой проблемы разработаны алгоритмы вывода траекторий (Trajectory Inference). Их задача — найти геометрический «скелет» данных в пространстве сниженной размерности (обычно после PCA, t-SNE или UMAP) и спроецировать клетки на этот скелет.

Подход на основе минимального остовного дерева (Slingshot)

Slingshot — алгоритм, работающий поверх уже существующей кластеризации и снижения размерности. Его архитектура включает два последовательных этапа:

Построение глобальной топологии. Вычисляются центроиды (центры масс) для каждого кластера в пространстве сниженной размерности. Между центроидами строится минимальное остовное дерево (Minimum Spanning Tree, MST). MST — это граф, соединяющий все узлы так, чтобы сумма длин ребер была минимальной, не образуя циклов. На этом этапе формируется грубый каркас: определяется, что кластер А переходит в кластер Б, а из кластера Б пути расходятся в В и Г.
Сглаживание через главные кривые. Биологические процессы не описываются ломаными линиями. Slingshot использует метод главных кривых (Principal Curves) — нелинейное обобщение метода главных компонент (PCA). Главная кривая — это гладкая линия, проходящая через «середину» облака точек так, чтобы каждая точка данных проецировалась на нее кратчайшим путем. Алгоритм итеративно изгибает отрезки MST, подгоняя их под локальную плотность клеток.

Качество работы Slingshot критически зависит от качества предварительного снижения размерности. Если алгоритм UMAP «разорвал» непрерывный континуум клеток на два изолированных острова из-за неверно подобранных гиперпараметров (например, слишком малого значения n_neighbors), Slingshot не сможет построить между ними достоверную траекторию.

Подход на основе графов (Monocle 3)

Семейство методов, к которому принадлежит Monocle 3, использует концепцию обратного вложения графов (Reversed Graph Embedding). Вместо того чтобы опираться на предварительную кластеризацию, алгоритм пытается напрямую выучить структуру главных графов (Principal Graphs) в пространстве UMAP.

Метод распределяет по облаку клеток узлы графа и итеративно оптимизирует их положение и связи между ними, минимизируя расстояние от каждой клетки до ближайшего ребра графа. Этот подход эффективнее справляется со сложными топологиями, включающими множественные ветвления и замкнутые циклы, что делает его незаменимым при анализе клеточного цикла или сложных иерархий гемопоэза.

Проблема корня (Root selection)

Фундаментальное ограничение большинства топологических методов: ни Slingshot, ни Monocle изначально не знают направления стрелы времени. Математически траектория от стволовой клетки к нейрону идентична траектории от нейрона к стволовой клетке. Алгоритм прокладывает маршрут, но не векторы движения.

Чтобы задать направление, исследователь обязан вручную указать «корень» (root) траектории. Это требует априорных биологических знаний. Выбирается кластер с максимальной экспрессией маркеров недифференцированных клеток (например, ген CD34 для гемопоэтических стволовых клеток), и этой точке присваивается значение псевдовремени, равное нулю.

Динамика экспрессии: поиск драйверов дифференцировки

После вычисления координаты псевдовремени фокус анализа смещается на биологическую интерпретацию: какие гены управляют движением клетки по найденной траектории?

Стандартные методы дифференциальной экспрессии из bulk RNA-Seq (DESeq2, edgeR) или базовые статистические тесты (U-критерий Манна-Уитни) здесь неприменимы. Они сравнивают дискретные группы (больные против здоровых, кластер А против кластера Б). В траекторном анализе переменная времени непрерывна. Необходимо найти гены, чья экспрессия статистически значимо меняется вдоль псевдовремени.

Для этой задачи применяются обобщенные аддитивные модели (Generalized Additive Models, GAM). В отличие от линейной регрессии, предполагающей монотонный рост или падение, GAM использует сглаживающие сплайны — кусочно-полиномиальные функции.

Базовое уравнение GAM в контексте псевдовремени выглядит так:

g(\mu_i) = \beta_0 + f(t_i)

Где $g$ — функция связи (обычно логарифмическая), $\mu_i$ — ожидаемая экспрессия гена в клетке $i$ , $\beta_0$ — базовый уровень экспрессии (интерсепт), а $f(t_i)$ — нелинейная гладкая функция (сплайн) от псевдовремени $t_i$ .

Модель проверяет нулевую гипотезу: функция $f(t_i)$ равна нулю на всем протяжении псевдовремени (экспрессия константна). Если ген значимо отклоняется от константы (например, сначала растет, достигает пика в середине траектории, а затем падает), он признается динамически экспрессирующимся.

Классический пример — бифуркация при гемопоэзе, когда мультипотентный предшественник дифференцируется либо в эритроцит, либо в клетку миелоидного ряда. Построив графики экспрессии вдоль двух расходящихся ветвей псевдовремени, можно наблюдать работу генных регуляторных сетей. На ранних этапах (до ветвления) экспрессируются оба транскрипционных фактора: GATA1 (драйвер эритропоэза) и PU.1 (драйвер миелопоэза). По мере продвижения к точке бифуркации начинается взаимный антагонизм. На ветви эритроцитов экспрессия GATA1 резко возрастает, подавляя PU.1, а на миелоидной ветви картина зеркально меняется. GAM-модели позволяют математически точно локализовать точку в псевдовремени, где происходит этот перелом.

RNA Velocity: векторное поле клеточных состояний

Топологические методы (Slingshot, Monocle) имеют концептуальный изъян: они опираются исключительно на транскриптомное сходство клеток. Эвристика «если клетки похожи, значит, они переходят друг в друга» работает не всегда. Революция в траекторном анализе произошла с появлением концепции RNA Velocity (РНК-скорость), извлекающей информацию о направлении развития из физики транскрипции.

Генная экспрессия — процесс с собственной кинетикой. РНК-полимераза синтезирует пре-мРНК, содержащую экзоны и интроны. Сплайсосома вырезает интроны, формируя зрелую мРНК, которая в дальнейшем транслируется и в конечном итоге деградирует.

Стандартные пайплайны (Cell Ranger) подсчитывают риды, картирующиеся на экзоны. Однако в сырых данных FASTQ всегда присутствует фракция (10–20%) интронных прочтений — молекул, захваченных до завершения сплайсинга.

Инструменты RNA Velocity (например, velocyto или scVelo) разделяют транскриптом клетки на две матрицы:

$U$ (Unspliced) — несплайсированные транскрипты (недавно синтезированные, содержат интроны).
$S$ (Spliced) — сплайсированные транскрипты (зрелые молекулы, только экзоны).

Кинетика описывается системой обыкновенных дифференциальных уравнений:

\frac{dU}{dt} = \alpha - \beta U

\frac{dS}{dt} = \beta U - \gamma S

Где $\alpha$ — скорость транскрипции, $\beta$ — скорость сплайсинга, $\gamma$ — скорость деградации зрелой мРНК.

Ключевой инсайт заключается в анализе фазового портрета (отношения $U$ к $S$ ) для каждого гена. В состоянии равновесия (steady state), когда экспрессия стабильна, синтез компенсирует распад: $\beta U = \gamma S$ . На графике, где по оси абсцисс отложено $S$ , а по оси ординат $U$ , равновесные клетки выстраиваются вдоль диагональной линии.

Если ген Neurog2 (маркер нейрогенеза) только начал активно транскрибироваться, баланс нарушается. Количество несплайсированной РНК ( $U$ ) резко возрастает, а зрелая РНК ( $S$ ) еще не накопилась. Точка клетки на фазовом портрете отклоняется выше равновесной диагонали. При выключении гена транскрипция падает ( $\alpha$ стремится к нулю), $U$ быстро исчезает за счет сплайсинга, а $S$ медленно деградирует. Точка отклоняется ниже диагонали.

Решая эти уравнения, алгоритм вычисляет вектор скорости для каждого гена в каждой клетке. Положительная скорость означает, что экспрессия гена в ближайшие часы вырастет, отрицательная — упадет. Сложив векторы всех генов, формируется многомерный вектор, указывающий будущее транскриптомное состояние клетки. Проекция этих векторов на UMAP создает векторное поле, где стрелки показывают биологические потоки. Важнейшее преимущество: RNA Velocity не требует ручного указания «корня» — направление выводится из внутренних кинетических законов.

Интеграция Bulk и Single-Cell данных в траекторном анализе

Важным этапом комплексного транскриптомного исследования является валидация scRNA-Seq траекторий независимыми методами. Данные bulk RNA-Seq, полученные из отсортированных (FACS) популяций клеток на известных стадиях развития, служат мощным инструментом верификации.

Интеграция происходит путем проецирования профилей bulk-секвенирования на пространство сниженной размерности scRNA-Seq (например, с помощью взаимных ближайших соседей — MNN, или канонического корреляционного анализа — CCA в Seurat). Если математически выведенная траектория верна, то bulk-образцы, отсортированные по маркерам ранних стадий, спроецируются в начало псевдовремени, а образцы поздних стадий — в конец. Такое кросс-платформенное подтверждение защищает исследование от артефактов вычислительных алгоритмов.

Ограничения и ловушки

Траекторный анализ требует жесткого критического контроля. Известная в биоинформатике шутка гласит: «Если скормить алгоритму мешок картошки, он найдет в нем траекторию дифференцировки».

Разрывы плотности. Алгоритмы предполагают непрерывность. Если промежуточные состояния быстротечны или выживаемость клеток на этапе перехода низка, алгоритм может провести ложный мост через пустоту между двумя биологически не связанными типами клеток.
Дискретные переходы. Клеточное слияние, реакция на острый стресс или апоптоз происходят скачкообразно. Натягивание непрерывной траектории на такие события искажает биологический смысл.
Батч-эффекты (Batch effects). Если клетки ранних стадий секвенированы в одной партии (batch), а поздних — в другой, технические различия могут сформировать ложную траекторию, отражающую не биологию, а разницу в протоколах пробоподготовки.
Нарушение кинетики в RNA Velocity. Модель предполагает постоянство скоростей сплайсинга ( $\beta$ ) и деградации ( $\gamma$ ). Однако вирусы или токсины могут целенаправленно блокировать сплайсосому. В таких случаях изменение отношения $U/S$ отражает глобальное нарушение метаболизма РНК, а не смену клеточного состояния, генерируя ложные векторные поля.

Интеграция мультиомиксных данных: совместный анализ и деконволюция bulk и single-cell RNA-Seq

Анализ транскриптома солидной опухоли из когорты в 500 пациентов показывает резкое, статистически значимое повышение экспрессии генов главного комплекса гистосовместимости (MHC-II) и интерферона-гамма после применения экспериментальной терапии. Делает ли препарат сами раковые клетки более иммуногенными, заставляя их синтезировать эти белки? Или же терапия просто спровоцировала массовую инфильтрацию опухоли Т-лимфоцитами и макрофагами, которые и принесли с собой эту РНК? Классический bulk RNA-Seq, усредняющий сигнал от миллионов клеток в гомогенизированном куске ткани, принципиально не способен ответить на этот вопрос. Сигнал об изменении внутриклеточной экспрессии и сигнал об изменении клеточного состава слиты воедино.

Современный транскриптомный анализ решает эту проблему через интеграцию двух парадигм. Single-cell RNA-Seq (scRNA-Seq) дает высочайшее клеточное разрешение, позволяя изучать гетерогенность ткани, но остается дорогим, технически зашумленным и часто ограничивается малыми выборками (например, 3 пациента против 3 контролей). Bulk RNA-Seq дешев, обладает огромной статистической мощностью на больших когортах и не страдает от артефактов диссоциации ткани, при которых хрупкие клетки (например, нейроны или адипоциты) разрушаются до попадания в секвенатор. Совместный анализ позволяет использовать scRNA-Seq как микроскоп для создания точного «словаря» клеточных типов, а bulk RNA-Seq — как телескоп для поиска статистически достоверных закономерностей на масштабе популяций.

Деконволюция: математическое разделение клеточных смесей

Процесс оценки пропорций различных типов клеток в смешанном (bulk) образце на основе референсных профилей экспрессии называется транскриптомной деконволюцией.

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

$B = S \times P$

Где:

$B$ (Bulk) — вектор или матрица экспрессии генов в смешанных образцах (известная величина, получаемая из сырых чтений bulk RNA-Seq после нормализации).
$S$ (Signature) — сигнатурная матрица, где столбцы — это типы клеток, строки — маркерные гены, а ячейки — средний уровень экспрессии гена в данном типе клеток (известная величина, извлекаемая из аннотированных данных scRNA-Seq).
$P$ (Proportions) — вектор пропорций клеточных типов в смешанном образце (неизвестная искомая величина, отражающая клеточный состав).

Поскольку матрица $B$ и матрица $S$ нам известны, задача сводится к нахождению вектора $P$ . В классической линейной алгебре это решается методом наименьших квадратов (Ordinary Least Squares, OLS). Однако в биологических системах применение чистого OLS часто выдает абсурдные результаты: алгоритм может решить, что наилучшее математическое приближение достигается, если доля В-клеток составит $40\%$ , а доля Т-клеток составит $-15\%$ .

Отрицательных количеств клеток в ткани не существует. Поэтому фундаментальной основой алгоритмов деконволюции является NNLS (Non-Negative Least Squares) — неотрицательный метод наименьших квадратов.

Целевая функция базового алгоритма NNLS выглядит так:

$\min ||B - S \times P||_2^2 \quad \text{при условии} \quad P_i \geq 0 \quad \text{и} \quad \sum P_i = 1$

Где:

$\min ||B - S \times P||_2^2$ — требование минимизировать квадрат евклидова расстояния (ошибку реконструкции) между реальным bulk-профилем и искусственным профилем, реконструированным из перемножения сигнатур на пропорции.
$P_i \geq 0$ — жесткое биологическое ограничение: доля любого $i$ -го клеточного типа должна быть больше или равна нулю.
$\sum P_i = 1$ — сумма долей всех клеточных типов должна равняться единице (или $100\%$ ), так как мы описываем замкнутую систему ткани.

Подготовка сигнатурной матрицы в Seurat и Scanpy

Точность деконволюции критически зависит от качества сигнатурной матрицы $S$ . Кастомизация пайплайна начинается именно здесь. Если включить в матрицу все 20 000 белок-кодирующих генов, алгоритм утонет в техническом шуме. Если включить гены, которые одинаково высоко экспрессируются в CD4+ и CD8+ Т-клетках, матрица станет коллинеарной — математически эти столбцы будут почти неразличимы, и алгоритм NNLS начнет случайным образом перебрасывать «веса» между ними, выдавая нестабильные пропорции.

Для создания надежной матрицы $S$ биоинформатики используют инструменты вроде Seurat (R) или Scanpy (Python). Стандартный путь включает:

Фильтрацию артефактов: удаление клеток с аномально высокой долей митохондриальных генов (свидетельство разрушения клетки) и рибосомальных генов. Пороги кастомизируются под ткань: для кардиомиоцитов $20\%$ митохондриальной РНК — норма, для лимфоцитов — признак апоптоза.
Кластеризацию и аннотацию: использование алгоритмов снижения размерности (UMAP/t-SNE) и графовой кластеризации (Leiden/Louvain) для выделения биологически осмысленных популяций.
Поиск маркерных генов: применение статистических тестов (например, Wilcoxon rank-sum test) для поиска генов, уникально экспрессирующихся в каждом кластере. Отбираются гены с высоким логарифмическим изменением кратности (Log2FC $> 1$ ) и строгим порогом скорректированного p-value.

Только эти высокоспецифичные маркерные гены формируют строки сигнатурной матрицы $S$ .

Продвинутые алгоритмы: CIBERSORTx и MuSiC

Базовый NNLS работает идеально только в симуляциях. В реальных данных возникают кросс-платформенные батч-эффекты (разница между технологиями секвенирования) и высокая межиндивидуальная вариабельность. Два наиболее востребованных инструмента решают эти проблемы разными путями.

CIBERSORTx: устойчивость к батч-эффектам и неизвестным клеткам

Оригинальный алгоритм CIBERSORT использовал метод опорных векторов для регрессии (Support Vector Regression, SVR) вместо обычного NNLS. SVR позволяет штрафовать модель не за все ошибки подряд, а только за те, которые выходят за рамки заданного порога толерантности, что делает алгоритм устойчивым к шуму в данных экспрессии.

Главное нововведение обновленной версии CIBERSORTx — способность генерировать кастомные сигнатурные матрицы непосредственно из пользовательских scRNA-Seq данных с автоматической коррекцией кросс-платформенного батч-эффекта (B-mode и S-mode коррекция). Алгоритм математически моделирует тот факт, что профиль экспрессии макрофага, захваченного в каплю микрофлюидной системы 10x Genomics, и макрофага, секвенированного в составе цельного куска ткани по протоколу Illumina TruSeq, систематически различаются из-за разной эффективности захвата мРНК.

Кроме того, CIBERSORTx элегантно решает проблему «неизвестного содержимого». В опухолевой ткани раковые клетки обладают уникальным, мутировавшим транскриптомом у каждого пациента. Создать универсальную сигнатуру раковой клетки невозможно. CIBERSORTx вводит концепцию «неизвестной фракции» (unassigned fraction). Алгоритм оценивает долю известных иммунных и стромальных клеток, не пытаясь принудительно распределить $100\%$ сигнала bulk-образца только по референсным сигнатурам, оставляя «остаточный» сигнал на долю опухоли.

MuSiC: учет кросс-субъектной вариативности

Алгоритм MuSiC (Multi-subject Single-cell deconvolution) обращает внимание на структуру самих scRNA-Seq данных. Обычно референсный датасет состоит из клеток, полученных от нескольких доноров, а не от одного идеального субъекта.

Рассмотрим пример построения сигнатуры бета-клеток поджелудочной железы. Ген INS (кодирующий инсулин) экспрессируется очень высоко, но у Донора 1 его средний уровень в бета-клетках составляет 10 000 каунтов, а у Донора 2 — 50 000 каунтов. Ген PDX1 экспрессируется на скромном уровне 500 каунтов, но этот уровень стабилен у всех доноров. Если использовать простое усреднение всех клеток (как делает большинство алгоритмов), ген INS будет доминировать в сигнатурной матрице, но его огромная межиндивидуальная дисперсия сделает деконволюцию нестабильной при применении к новым пациентам.

MuSiC использует взвешенный NNLS (Weighted NNLS). Алгоритм придает больший статистический вес тем маркерным генам, которые демонстрируют низкую дисперсию между биологическими субъектами в scRNA-Seq датасете. Гены, которые сильно «скачут» от пациента к пациенту, пенализируются (их вес снижается), даже если они являются классическими маркерами типа клеток с высокой экспрессией. Это делает перенос вычисленных пропорций на независимую bulk-когорту значительно более точным и биологически релевантным.

Ловушка псевдорепликации и Pseudobulk анализ

Если деконволюция — это проецирование знаний от единичных клеток на bulk-данные, то существует и обратный процесс: применение строгих статистических стандартов bulk-анализа к scRNA-Seq данным. Эта концепция называется псевдобалк (pseudobulk) агрегацией.

Рассмотрим классический дизайн эксперимента: анализируются 3 мыши дикого типа (WT) и 3 мыши с нокаутом гена (KO). Из каждой мыши с помощью scRNA-Seq выделено и отсеквенировано по 5 000 фибробластов. Итого получено 30 000 клеток. Задача — найти гены, дифференциально экспрессирующиеся в фибробластах при нокауте целевого гена.

Если применить стандартный клеточный тест, встроенный в Seurat (например, U-критерий Манна-Уитни), ко всем клеткам напрямую, сравнивая 15 000 клеток WT против 15 000 клеток KO, исследователь совершит грубейшую статистическую ошибку — псевдорепликацию.

Алгоритм будет исходить из того, что размер выборки $N = 30 000$ независимых наблюдений. При таком гигантском $N$ любое, даже самое крошечное техническое отклонение (например, одна мышь из группы WT случайно получила чуть больше стресса при диссоциации ткани, что слегка повысило экспрессию генов теплового шока) будет признано статистически значимым с $p < 10^{-100}$ . Однако клетки, извлеченные из одной мыши, не являются независимыми биологическими наблюдениями — они делят общий генетический фон и общую среду. Реальный размер биологической выборки в этом эксперименте — $N = 6$ (3 против 3).

Механика псевдобалк агрегации и работа с DESeq2

Чтобы получить истинные биологические выводы и избежать лавины ложноположительных результатов, данные необходимо агрегировать на уровне биологических реплик.

Процесс кастомизации пайплайна под pseudobulk состоит из следующих шагов:

Клетки кластеризуются и аннотируются стандартными методами в Seurat/Scanpy.
Внутри каждого биологического образца изолируются все клетки конкретного типа (например, выделяются только фибробласты Мыши 1, затем фибробласты Мыши 2 и т.д.).
Сырые чтения (raw counts) каждого гена суммируются по всем клеткам данного типа внутри каждого образца.
В результате для фибробластов формируется классическая матрица каунтов, где строки — гены, а столбцы — биологические образцы (WT_1, WT_2, WT_3, KO_1, KO_2, KO_3).
К этой матрице применяются классические, проверенные временем инструменты для bulk RNA-Seq, такие как DESeq2 или edgeR.

Суммирование сырых каунтов (а не усреднение уже нормализованных значений) критически важно. Алгоритм DESeq2 ожидает на вход строго целочисленные данные, так как использует отрицательное биномиальное распределение для моделирования связи между средним значением экспрессии и дисперсией.

Математически модель DESeq2 для каунтов гена $i$ в образце $j$ описывается как:

$K_{ij} \sim NB(\mu_{ij}, \alpha_i)$

Где:

$K_{ij}$ — наблюдаемое количество сырых чтений (counts).
$NB$ — отрицательное биномиальное распределение (Negative Binomial), идеально описывающее передисперсию (overdispersion) в данных секвенирования.
$\mu_{ij}$ — ожидаемое среднее значение, которое масштабируется на size factor (фактор размера библиотеки образца $j$ ).
$\alpha_i$ — параметр дисперсии для гена $i$ , оцениваемый на основе биологических реплик.

Суммируя каунты, мы сохраняем информацию о глубине секвенирования каждого «псевдообразца». Если в WT_1 было 5000 фибробластов, а в WT_2 — только 1000, суммарный профиль WT_1 будет иметь в 5 раз больше каунтов. DESeq2 автоматически учтет это через вычисление size factors, придавая более глубоко отсеквенированным образцам правильный статистический вес, в то время как простое усреднение уничтожило бы эту информацию.

Граничные случаи и фильтрация редких популяций

Псевдобалк анализ требует особой осторожности при работе с редкими клеточными популяциями. Если в образце WT_1 захвачено 2000 макрофагов, а в образце WT_2 — всего 5 макрофагов (из-за проблем с выживаемостью клеток при ферментативной обработке ткани), суммированный профиль WT_2 будет крайне разреженным, состоящим в основном из нулей.

Кастомизация пайплайна здесь заключается во введении жестких порогов фильтрации:

Исключать из анализа образцы, в которых представлено менее заданного числа клеток целевого типа (обычно порог составляет от 10 до 30 клеток).
Если после исключения некачественных образцов в одной из экспериментальных групп остается менее 2-3 биологических реплик, дифференциальный анализ для данного типа клеток проводить нельзя. Алгоритму DESeq2 не хватит степеней свободы для достоверной оценки внутригрупповой дисперсии $\alpha_i$ , что приведет к некорректным p-values.

Синтез подходов в реальном исследовании

Интеграция мультиомиксных данных на практике часто представляет собой замкнутый аналитический цикл, где методы дополняют друг друга.

На первом этапе исследователь использует крупную когорту bulk RNA-Seq (например, 1000 образцов транскриптомов опухолей из базы TCGA), чтобы найти модули генов, достоверно коррелирующие с выживаемостью пациентов. На втором этапе проводится scRNA-Seq на малой, но детально охарактеризованной выборке (5-10 пациентов), чтобы составить атлас микроокружения опухоли и понять, какие именно клеточные субпопуляции физически экспрессируют эти гены выживаемости. На третьем этапе (деконволюция) сигнатурные матрицы, извлеченные из scRNA-Seq, применяются к изначальной тысяче bulk-образцов. Это позволяет доказать, что выживаемость коррелирует не просто с абстрактным списком генов, а с физической инфильтрацией опухоли конкретным подтипом клеток (например, истощенными CD8+ Т-клетками). На четвертом этапе (псевдобалк) внутри scRNA-Seq датасета проводится строгий дифференциальный анализ между пациентами, ответившими и не ответившими на терапию, чтобы найти новые терапевтические мишени непосредственно внутри выявленной ключевой клеточной популяции.

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

Обеспечение воспроизводимости исследований и кастомизация биоинформатических пайплайнов под научные задачи

Вы открываете свой же Jupyter-ноутбук или R-скрипт полугодовой давности, чтобы добавить в анализ один контрольный образец по просьбе рецензента журнала. Нажимаете «Run All» и ожидаете увидеть те же самые UMAP-кластеры и списки дифференциально экспрессирующихся генов, лишь слегка скорректированные новыми данными. Вместо этого скрипт падает на третьем шаге: функция из пакета Seurat не может прочитать объект, так как за это время библиотека обновилась с мажорной версии 3 на 4, полностью изменив внутреннюю структуру классов. Когда вы чудом откатываете версию и чините ошибку, итоговое количество кластеров клеток меняется с восьми на одиннадцать. Биологические выводы статьи рушатся.

Эта ситуация — классическое проявление кризиса воспроизводимости в вычислительной биологии. Итоговый результат транскриптомного анализа зависит не только от исходных FASTQ-файлов, но и от сотен скрытых переменных: версий операционной системы, системных библиотек, генераторов псевдослучайных чисел и порядка выполнения команд. Переход от интерактивного запуска команд к созданию надежных, воспроизводимых и масштабируемых конвейеров (пайплайнов) — это граница, отделяющая базовое владение биоинформатикой от профессионального инжиниринга данных.

Иллюзия работающего скрипта и уровни изоляции

Традиционный подход начинающего исследователя заключается в написании длинного bash-скрипта (например, run_analysis.sh), в котором последовательно вызываются инструменты: контроль качества, тримминг, выравнивание, подсчет каунтов. Главная уязвимость такого подхода — жесткая привязка к окружению, в котором этот скрипт был написан.

Инструменты для RNA-Seq зависят от специфических версий компиляторов (C/C++), библиотек линейной алгебры и интерпретаторов (Python, R). Обновление системной библиотеки libcurl или zlib на сервере может незаметно изменить поведение алгоритма или сломать компиляцию пакета DESeq2. Чтобы изолировать анализ от нестабильности внешней среды, применяются два уровня управления зависимостями.

Пакетные менеджеры: Conda и Mamba

Первый уровень изоляции — создание виртуальных окружений. Менеджер Conda позволяет описать все необходимые программы в одном текстовом файле environment.yml.

name: rna_seq_env
channels:
  - conda-forge
  - bioconda
dependencies:
  - fastqc=0.11.9
  - star=2.7.10a
  - subread=2.0.3
  - r-deseq2=1.38.3

При развертывании этого окружения Conda скачивает бинарные файлы программ и их зависимости в изолированную директорию, не затрагивая системные файлы. Однако базовая Conda использует классический SAT-солвер (алгоритм решения задачи выполнимости булевых формул) для разрешения графа зависимостей. При установке сложных биоинформатических сред с сотнями пакетов этот процесс может занимать часы.

Mamba — это drop-in замена для Conda, написанная на C++. Она использует более современные алгоритмы разрешения зависимостей (на базе библиотеки libsolv), сокращая время установки сложных сред с часов до минут, при этом сохраняя полную совместимость с форматом environment.yml.

Тем не менее, виртуальные окружения не решают проблему базовой операционной системы. Если пакет был скомпилирован под Ubuntu 18.04, он может вести себя иначе или не запуститься на CentOS 7 из-за разницы в базовых системных вызовах ядра Linux (glibc).

Контейнеризация: Docker и Singularity

Второй, абсолютный уровень изоляции — контейнеризация. Контейнер упаковывает в единый образ не только биоинформатические программы, но и всю файловую систему операционной системы (кроме самого ядра).

Docker стал стандартом де-факто в IT-индустрии. Процесс создания образа описывается в Dockerfile, где фиксируется базовая ОС (например, Debian 11), установка системных библиотек и биоинформатического софта. Образ собирается один раз, получает уникальный криптографический хэш (SHA256) и может быть запущен на любом компьютере мира с гарантией побитового совпадения среды.

В биоинформатике возникает фундаментальный нюанс: архитектура Docker требует прав суперпользователя (root) для запуска фонового демона. На высокопроизводительных вычислительных кластерах (HPC), где анализируются терабайты scRNA-Seq данных, администраторы никогда не дадут пользователям root-права. Если пользователь имеет доступ к демону Docker, он может примонтировать корневую файловую систему хоста внутрь контейнера и получить полный несанкционированный контроль над суперкомпьютером.

Решением стала технология Singularity (ныне развивающаяся под именем Apptainer). Singularity использует пространства имен пользователя (user namespaces) ядра Linux, что позволяет запускать контейнеры от имени обычного пользователя без повышения привилегий. Система прозрачно монтирует домашние директории кластера внутрь контейнера. Исследователь может создать образ в Docker на локальном ноутбуке, конвертировать его в формат .sif (Singularity Image Format) и безопасно запустить на кластере института.

Управление потоками данных: от bash к Nextflow

Даже если все программы упакованы в контейнеры, монолитный bash-скрипт остается уязвимым к сбоям. При обработке 100 образцов scRNA-Seq выравнивание 99 образцов может пройти успешно, а на сотом сервер перезагрузится из-за нехватки оперативной памяти. При использовании bash-скрипта потребуется вручную комментировать строки с успешными образцами, чтобы перезапустить только упавший. В сложных конвейерах ручное вмешательство неизбежно приводит к потере данных или путанице.

Современная транскриптомика опирается на системы управления рабочими процессами (Workflow Managers), лидерами среди которых являются Nextflow (на базе языка Groovy) и Snakemake (на базе Python).

Направленный ациклический граф (DAG)

В основе систем управления лежит математическая концепция направленного ациклического графа. Граф можно выразить формулой $G = (V, E)$ , где $G$ — это сам граф вычислений, $V$ (вершины) — вычислительные процессы (например, выравнивание STAR, подсчет каунтов FeatureCounts), а $E$ (ребра) — потоки данных (например, FASTQ- или BAM-файлы), передаваемые от одного процесса к другому. Ацикличность означает, что данные текут строго вперед от сырых чтений к финальной матрице, без бесконечных циклов.

Nextflow автоматически строит этот граф на основе написанного кода. Разработчику не нужно указывать, какие задачи запускать параллельно. Достаточно определить входы и выходы процессов. Если на вход подается 100 пар FASTQ-файлов, Nextflow анализирует граф $G$ , понимает, что процесс FASTQC для первого образца не зависит от второго, и автоматически запускает 100 параллельных задач (jobs), отправляя их в планировщик кластера (например, SLURM или PBS).

Кэширование и реентерабельность

Главная вычислительная сила Nextflow кроется в механизме возобновления (resume). Для каждой задачи (вершины $V$ ) Nextflow вычисляет уникальный 128-битный хэш, который зависит от трех компонентов:

Содержимого и метаданных входных файлов (ребер $E$ ).
Текста самого скрипта внутри процесса.
Версии используемого контейнера или окружения.

Если задача завершилась успешно, ее результат сохраняется в скрытой рабочей директории. Если конвейер упал на этапе сборки итоговой матрицы экспрессии, достаточно запустить ту же команду с флагом -resume. Nextflow заново вычисляет хэши графа: выравнивание и тримминг не изменились, поэтому система мгновенно берет их результаты из кэша (создавая символические ссылки) и продолжает вычисления ровно с упавшего узла. Это экономит тысячи часов машинного времени и бюджет лаборатории.

Экосистема nf-core: стандартизация аналитики

Исторически каждая исследовательская группа писала свои конвейеры. Сравнение двух статей по scRNA-Seq часто превращалось в сравнение артефактов самописных скриптов, а не реальной биологии. Разные параметры тримминга или версии аннотаций генома давали расхождения в списках дифференциально экспрессирующихся генов до 30%.

Проект nf-core изменил эту парадигму. Это глобальное сообщество биоинформатиков, разрабатывающее стандартизированные, строго рецензируемые пайплайны на Nextflow. Пайплайны nf-core/rnaseq (для bulk) и nf-core/scrnaseq (для единичных клеток) признаны золотым стандартом индустрии.

Преимущества использования nf-core:

Встроенная контейнеризация: Каждый шаг пайплайна автоматически подтягивает нужный, протестированный Docker/Singularity контейнер. Отпадает необходимость локальной установки софта.
Агностичность к инфраструктуре: Пайплайн идентично работает на локальном сервере (с флагом -profile docker), на университетском кластере (с флагом -profile singularity) и на облачных платформах AWS/Google Cloud (с флагом -profile awsbatch).
Best Practices: В nf-core/rnaseq уже встроены контроль качества (FastQC), удаление адаптеров (TrimGalore), псевдовыравнивание (Salmon), классическое выравнивание (STAR), оценка контаминации рибосомальной РНК и автоматическая генерация сводного интерактивного отчета MultiQC.

Запуск полного транскриптомного анализа сводится к одной команде, где передается таблица с путями к сырым данным и ссылки на референсы:

nextflow run nf-core/rnaseq \
  -profile singularity \
  --input samplesheet.csv \
  --fasta genome.fa \
  --gtf annotation.gtf \
  --outdir results

Кастомизация пайплайнов под научные задачи

Несмотря на мощь nf-core, реальные научные задачи часто требуют выхода за рамки стандартов. Исследование вирусной инфекции на уровне единичных клеток — типичный пример. Стандартный пайплайн nf-core/scrnaseq выровняет риды на геном человека, но проигнорирует вирусные транскрипты, выбросив их как некартированные (unmapped). В результате теряется информация о том, какие именно клетки инфицированы.

Кастомизация требует понимания того, как перехватить поток данных или изменить конфигурацию процесса, не разрушая базовый код пайплайна.

Уровень 1: Модификация конфигурации

Самый безопасный способ кастомизации — использование пользовательских конфигурационных файлов (custom.config). В случае с вирусом нет необходимости переписывать код nf-core. Достаточно создать объединенный химерный референсный геном (FASTA человека + FASTA вируса) и объединенную аннотацию (GTF).

Если вирусные транскрипты имеют нестандартную структуру (например, перекрывающиеся рамки считывания), потребуется изменить параметры выравнивателя STARsolo, который работает внутри пайплайна. В custom.config переопределяются аргументы конкретного процесса (вершины графа):

process {
    withName: 'STARSOLO' {
        ext.args = '--clip3pAdapterSeq AAAAAA --outFilterMatchNmin 30'
    }
}

Передав этот файл при запуске (-c custom.config), исследователь инжектирует специфические параметры в стандартный конвейер. Сам конвейер остается нетронутым, что позволяет легко обновлять его до новых версий в будущем.

Уровень 2: Инъекция собственных модулей

Более глубокий сценарий кастомизации возникает при работе со сложными тканями. При секвенировании единичных клеток (например, на платформе 10x Genomics) часть клеток разрушается в процессе пробоподготовки. Их РНК выливается в общий раствор и попадает в капли с целыми клетками. Это фоновая РНК (ambient RNA). Если ткань богата РНКазами или высокоэкспрессируемыми генами (например, ткань печени или поджелудочной железы), фоновый шум может создать ложные сигналы экспрессии во всех кластерах.

Стандартный пайплайн заканчивается на выдаче сырой матрицы каунтов. Чтобы очистить данные, матрицу необходимо пропустить через специализированный алгоритм (например, CellBender, который использует глубокие генеративные модели для вычитания фона) до того, как данные попадут во вторичный анализ (Seurat или Scanpy).

В Nextflow (с синтаксисом DSL2) пайплайны строятся из независимых модулей. Можно написать собственный скрипт-обертку, который импортирует стандартный пайплайн как подпрограмму, перехватывает его выходной канал (channel) с матрицами и направляет в кастомный процесс:

include { SCRNASEQ } from './workflows/scrnaseq'
include { CELLBENDER } from './modules/local/cellbender'

workflow {
    // 1. Запускаем стандартный пайплайн nf-core
    SCRNASEQ ()

    // 2. Перехватываем сырые матрицы и отправляем в очистку
    CELLBENDER ( SCRNASEQ.out.raw_matrices )
}

Такой модульный подход гарантирует, что базовая обработка (выравнивание, демультиплексирование) остается строго стандартизированной, а научная новизна (очистка от фона специфическим методом) элегантно наслаивается поверх.

Контроль случайности: Random Seeds

Важнейший аспект воспроизводимости, о котором часто забывают при переходе от первичной обработки (FASTQ $\rightarrow$ Матрица) к статистическому анализу (Матрица $\rightarrow$ Биологический смысл) — это стохастические алгоритмы.

Многие методы нелинейного снижения размерности (t-SNE, UMAP) и графовой кластеризации (Louvain, Leiden) опираются на генераторы псевдослучайных чисел. Например, UMAP использует стохастический градиентный спуск для оптимизации расположения точек в двумерном пространстве. Начальные координаты точек задаются случайным образом. Если запустить UMAP на одной и той же матрице экспрессии дважды без фиксации параметров, облака клеток будут выглядеть по-разному, а кластеры могут получить другие порядковые номера или слегка изменить границы.

Чтобы зафиксировать результат, необходимо явно задавать «зерно» (seed) генератора случайных чисел в скриптах (например, set.seed(42) в R или random_state=42 в Python-библиотеках). В профессиональных пайплайнах этот параметр выносится в глобальный конфигурационный файл, чтобы любой исследователь, запустивший код на другом континенте, получил математически идентичные координаты клеток на графике.

Версионирование кода и принцип FAIR

Воспроизводимость немыслима без контроля версий. Git используется для фиксации истории аналитического проекта. Однако критической ошибкой новичков является попытка коммитить в Git сами данные (матрицы экспрессии или BAM-файлы). Git предназначен исключительно для текстовых файлов (кода). Бинарные файлы раздувают репозиторий и делают его неработоспособным.

Современный стандарт публикации транскриптомных данных подчиняется принципам FAIR (Findable, Accessible, Interoperable, Reusable — Находимые, Доступные, Интероперабельные, Повторно используемые). Чтобы исследование считалось FAIR, необходимо выполнить три условия:

Сырые данные (FASTQ), обработанные матрицы каунтов и подробные метаданные образцов загружаются в специализированные публичные репозитории (GEO — Gene Expression Omnibus, или ArrayExpress).
Код анализа (Nextflow пайплайны, R/Python скрипты для построения графиков) публикуется в открытых репозиториях (GitHub, GitLab, Zenodo).
В репозитории с кодом фиксируются точные версии образов Singularity/Docker и указываются хэши коммитов используемых пайплайнов.

Фраза в статье «Анализ проводился с использованием Seurat v4» невоспроизводима, так как внутри мажорной версии v4 было множество минорных обновлений. Эталон воспроизводимой науки звучит так: «Анализ проведен с помощью пайплайна nf-core/scrnaseq v2.1.0 с использованием образа docker://nfcore/scrnaseq:2.1.0. Код вторичного анализа доступен в репозитории [Ссылка на GitHub] (commit 8f3a2b), seed=123. Сырые данные и матрицы доступны в GEO под номером GSEXXXXXX».

Фундаментальная база и рекомендуемые источники

Для глубокого погружения в тему воспроизводимости и кастомизации биоинформатических пайплайнов рекомендуется опираться на следующие проверенные материалы.

Научные статьи и стандарты:

Wilkinson, M. D., et al. (2016). "The FAIR Guiding Principles for scientific data management and stewardship." Scientific Data. — Оригинальная публикация, заложившая основы концепции FAIR. Обязательна к прочтению для понимания стандартов публикации данных.
Ewels, P. A., et al. (2020). "The nf-core framework for community-curated bioinformatics pipelines." Nature Biotechnology. — Статья создателей nf-core, подробно описывающая архитектуру и философию стандартизированных пайплайнов.

Учебники и официальные мануалы:

Nextflow Documentation (nextflow.io/docs) — Исчерпывающее руководство по синтаксису Groovy/DSL2, управлению каналами (channels) и настройке профилей для различных кластеров.
Single-cell Best Practices (sc-best-practices.org) — Современный, постоянно обновляемый онлайн-учебник от ведущих биоинформатиков (включая лабораторию Theis Lab). Содержит детальные разборы математики и кода для каждого этапа scRNA-Seq.
Antao, T. "Bioinformatics with Python Cookbook" — Отличная книга для понимания интеграции Python-инструментов с системными вызовами и контейнерами.

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

nf-core training (на YouTube и сайте nf-co.re) — Серия воркшопов и обучающих видео, где разработчики пошагово показывают запуск, кастомизацию и написание собственных модулей.
Курсы Института Биоинформатики (Bioinformatics Institute) — Русскоязычные курсы на платформах Stepik и YouTube. Особенно полезны модули по введению в Linux, bash-скриптингу и основам RNA-Seq.
Официальные виньетки Seurat (satijalab.org/seurat) и туториалы Scanpy (scanpy.readthedocs.io) — Базовые руководства по вторичному анализу, где наглядно показана важность контроля параметров (включая random seeds) при кластеризации и снижении размерности.

Внедрение контейнеров и менеджеров рабочих процессов требует начальных инвестиций времени в изучение синтаксиса и концепций. Однако этот барьер многократно окупается при первой же необходимости отмасштабировать анализ с одного образца на сотню, или когда через год потребуется заново сгенерировать графики, будучи абсолютно уверенным, что результаты не исказятся из-за незаметного обновления системной библиотеки.

Дизайн биологического эксперимента и технологии подготовки библиотек для RNA-Seq

Анатомия биологического эксперимента

Биологические и технические повторности

Эффект партии (Batch Effect)

Качество РНК: фундамент библиотеки

Стратегии обогащения: избавление от балласта

Poly-A селекция (Обогащение мРНК)

Ribo-Zero (Истощение рРНК)

Сравнительная таблица стратегий

Синтез кДНК и проблема направленности (Strandedness)

Параметры секвенирования: глубина и длина рида

Рекомендуемые ресурсы для углубленного изучения

Первичная обработка сырых данных: контроль качества и путь от FASTQ к матрице экспрессии

Первичная обработка сырых данных: контроль качества и путь от FASTQ к матрице экспрессии

Анатомия сырых данных: формат FASTQ

Контроль качества и «стрижка» ридов (Trimming)

Сплайсированное выравнивание: поиск координат на геноме

Форматы SAM/BAM и язык CIGAR

Альтернативный путь: псевдовыравнивание

Квантификация: от координат к матрице экспрессии

Рекомендуемые ресурсы для углубленного изучения

Статистический анализ bulk RNA-Seq: методы нормализации и оценка дифференциальной экспрессии генов

Две проблемы сырых каунтов: глубина и композиция

Почему TPM и RPKM не подходят для дифференциальной экспрессии

Алгоритм Median of Ratios (DESeq2)

Статистическая природа данных: Отрицательное биномиальное распределение

Проблема малых выборок и сжатие дисперсии (Shrinkage)

Множественное тестирование и контроль FDR

Log2 Fold Change и независимая фильтрация

Рекомендуемые ресурсы для углубленного изучения

Функциональная аннотация и интерпретация результатов: анализ обогащения GO, KEGG и GSEA

Функциональная аннотация и интерпретация результатов: анализ обогащения GO, KEGG и GSEA

Базы данных: формализация биологического знания

Gene Ontology (GO): строгая иерархия терминов

KEGG и Reactome: метаболические и сигнальные пути

MSigDB: золотой стандарт для GSEA

Over-Representation Analysis (ORA): поиск аномальных пересечений

Критическая ошибка ORA: неправильный выбор фона

Ограничения ORA

Gene Set Enrichment Analysis (GSEA): анализ без порогов

Шаг 1: Ранжирование генома

Шаг 2: Алгоритм Running Sum

Шаг 3: Оценка результатов (ES, NES и Leading Edge)

Сравнение подходов: что выбрать?

Борьба с избыточностью результатов

Рекомендуемые ресурсы для углубленного изучения

Основы single-cell RNA-Seq: технологические платформы, захват клеток и специфическая предобработка данных

Основы single-cell RNA-Seq: технологические платформы, захват клеток и специфическая предобработка данных

Технологические платформы: от лунки к капле

Планшетные методы (Smart-seq2 / Smart-seq3)

Капельная микрофлюидика (10x Genomics, Drop-seq)

Математика инкапсуляции и проблема дуплетов

Баркодирование: идентификация клеток и молекул

Cell Barcode (Клеточный баркод)

Unique Molecular Identifier (UMI)

Специфическая предобработка данных

Идентификация реальных клеток (Cell Calling)

Контроль качества на уровне клеток (QC)

Рекомендуемые ресурсы для углубленного изучения

Математическое снижение размерности и алгоритмы кластеризации в анализе единичных клеток

Отбор признаков: фильтрация шума и поиск вариабельности

Метод главных компонент (PCA) как фундамент анализа

Нелинейное снижение размерности: топология и визуализация

t-SNE и решение проблемы скученности

UMAP: алгебраическая топология и глобальная структура

Графовая кластеризация: от матриц к сетям

Шаг 1: Построение KNN-графа

Шаг 2: Переход к SNN-графу (Shared Nearest Neighbor)

Шаг 3: Оптимизация модулярности (Louvain и Leiden)

Управление детализацией: параметр Resolution

Материалы для углубленного изучения математического аппарата

Идентификация клеточных типов, поиск маркерных генов и автоматизированная аннотация популяций

Идентификация клеточных типов, поиск маркерных генов и автоматизированная аннотация популяций

Статистика поиска маркеров: почему bulk-методы здесь ломаются

Анатомия идеального маркерного гена

Ручная аннотация и паттерны визуализации

Автоматизированная аннотация: перенос ярлыков с референса

Анатомия «мусорных» и переходных кластеров

Рекомендуемые источники и материалы для изучения

Траекторный анализ и моделирование динамики генной экспрессии в процессах клеточной дифференцировки