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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

!Интерактивный симулятор электрофореграммы РНК: слайдер деградации плавно трансформирует профиль от идеальных пиков 28S/18S до горба коротких фрагментов, а показатель RIN меняется от 10.0 до 1.0 с цветовой индикацией качества.

Показатель RIN диктует, какую технологию подготовки библиотеки допустимо использовать. При доступны любые методы. Если 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 селекция требует .

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

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

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

Позволяет анализировать любые типы РНК, включая некодирующие и гистоновые мРНК.

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

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

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

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

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

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

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

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

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

!Схема подготовки библиотеки РНК-секвенирования с использованием метода направленного расщепления cДНК.

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

Синтез первой цепи: Обратная транскриптаза строит первую цепь кДНК на матрице РНК, используя стандартные дезоксирибонуклеотиды (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 прибор читает фрагмент с обоих концов (например, нуклеотидов).

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

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

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

Сборка транскриптома de novo (без референсного генома), где парные концы помогают алгоритмам (например, Trinity) строить длинные контиги.

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

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

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

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

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

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

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

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

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

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

| Значение | Вероятность ошибки () | Точность прочтения | | :--- | :--- | :--- | | 10 | 1 к 10 () | 90.0% | | 20 | 1 к 100 () | 99.0% | | 30 | 1 к 1000 () | 99.9% | | 40 | 1 к 10000 () | 99.99% |

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

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

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

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

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

Контаминация адаптерами. Библиотеки для секвенирования содержат синтетические адаптеры на концах фрагментов. Если исходный биологический фрагмент кДНК оказался короче, чем заданная длина чтения (например, фрагмент 100 пар оснований при режиме секвенирования 150 пар оснований), полимераза прочитает весь биологический фрагмент и продолжит синтез, читая адаптер на противоположном конце (read-through). Попытка найти эту химерную последовательность в референсном геноме обречена на провал.

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

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

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

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

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

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

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

!Анимация сплайсированного выравнивания рида на геном по алгоритму STAR MMP: рид прикладывается к Экзону 1, упирается в границу интрона, правая часть перелетает через интрон и выравнивается на Экзон 2, после чего появляется пунктирная дуга сплайсинга.

Алгоритм 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-строки при выравнивании рида на референсный геном с учетом интронов.

Ключевые операторы 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, а размер равен 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...), а на пересечении стоят целые числа: сколько уникальных фрагментов кДНК было приписано данному гену в данном образце.

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

3. Статистический анализ 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 ридов.

В опыте: Ген В заберет на себя подавляющее большинство ридов (около ). На гены А и Б останется суммарно около 34 ридов (по 17 на каждый).

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

!Искажение относительной экспрессии генов при нормализации данных в эксперименте bulk RNA-Seq.

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

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

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

Где:

— нормализованное значение для гена .

— количество сырых прочтений (counts), картированных на ген .

— эффективная длина гена в килобазах.

— сумма отношений каунтов к длине для всех генов в данном образце.

Суть TPM в том, что сумма всех значений TPM в любом образце всегда равна ровно (одному миллиону). Это позволяет корректно сравнивать гены внутри одного образца: если 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.

!Пошаговая визуализация алгоритма Median of Ratios: от исходной матрицы каунтов через псевдо-референс и отношения к нормализованным значениям, с наглядным показом устойчивости медианы к выбросам.

Почему используется именно медиана, а не среднее арифметическое? Медиана абсолютно нечувствительна к экстремальным выбросам. Возвращаясь к нашему примеру с токсином: если ген В вырос в 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-распределении описывается формулой:

Где:

— ожидаемая дисперсия каунтов для гена.

— среднее значение экспрессии гена.

— параметр дисперсии (dispersion parameter), отражающий биологическую вариабельность.

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

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

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

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

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

Сначала рассчитывается индивидуальная оценка дисперсии для каждого из 20 000 генов по отдельности.

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

Наконец, индивидуальные оценки дисперсии "подтягиваются" (shrink) к этой кривой.

!Снижение изменчивости дисперсии генов путем приближения к тренду в анализе RNA-Seq.

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

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

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

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

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

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

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

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

Где:

— скорректированное -value (оценка FDR).

— исходное -value гена.

— общее количество протестированных генов.

— ранг гена в отсортированном списке (от 1 — самый значимый, до — наименее значимый).

| Ген | Сырое -value | Ранг () | Всего генов () | Расчет | Итог | | :--- | :--- | :--- | :--- | :--- | :--- | | Ген X | 0.0001 | 1 | 10000 | | 1.0000 (ограничивается сверху) -> 0.0001 | | Ген Y | 0.0010 | 2 | 10000 | | 0.0050 | | Ген Z | 0.0500 | 100 | 10000 | | 5.0000 (ограничивается 1.0) |

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

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

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

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

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

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

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

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

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

Функциональная аннотация и интерпретация результатов: анализ обогащения 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). Термины организованы от широких понятий к узким, при этом один термин может иметь несколько «родителей».

!Структура направленного ациклического графа Gene Ontology

> Правило истинного пути (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). Его логика опирается на пересечение множеств: мы берем список генов, которые достоверно изменили экспрессию (например, прошли жесткие фильтры и ), и проверяем, нет ли среди них аномально высокой доли генов из какого-то известного биологического пути.

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

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

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

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

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

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

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

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

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

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

Зависимость от жестких порогов. Чтобы получить список DEGs, мы обязаны выбрать субъективные пороги. Ген с попадет в анализ, а ген с будет выброшен, хотя биологически между ними нет никакой разницы.

Игнорирование силы эффекта. ORA относится к списку DEGs как к плоскому бинарному множеству (ген либо есть в списке, либо нет). Ген, экспрессия которого изменилась в 100 раз, имеет ровно такой же вес при пересечении множеств, как ген, изменившийся в 1.5 раза.

Потеря слабых, но скоординированных сигналов. Если в метаболическом пути экспрессия всех 50 ферментов синхронно выросла на 15% (что недостаточно для преодоления порога FDR для каждого отдельного гена, но биологически приведет к мощному усилению синтеза метаболита), ORA вообще не увидит этот путь. Ни один ген не попадет в выборку, и путь будет пропущен.

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

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

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

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

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

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

Где возвращает для генов с растущей экспрессией и для генов с падающей. Функция превращает маленькие p-value в большие положительные числа. Например, если ген имеет и , его метрика . Он отправится в самый низ списка. В середине списка сгрудятся гены с (их метрика ) — это биологический шум, не отреагировавший на стимул.

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

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

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

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

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

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

!Интерактивная визуализация алгоритма Running Sum из GSEA: анимированная кривая обогащения строится в реальном времени, показывая разницу между активированным и случайным путём.

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

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

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

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

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

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

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

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

Независимо от того, используете вы 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 главных биологических тем, затронутых в вашем эксперименте.

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

Основы 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'-конец). Молекула мРНК захватывается за поли-А хвост, и секвенируется лишь короткий фрагмент рядом с ним, что делает невозможным полноценный анализ сплайсинга.

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

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

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

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

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

Пустые капли ():

Одиночные клетки ():

Две и более клеток ():

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

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

Пустые капли:

Одиночные клетки:

Дуплеты:

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

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

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

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

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

Unique Molecular Identifier (UMI)

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

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

!Использование UMI-тегов для идентификации и корректного учета ПЦР-дупликатов в данных scRNA-seq.

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

Длина UMI в 12 нуклеотидов дает уникальных комбинаций. Риск того, что две разные молекулы одного и того же гена в одной клетке случайно свяжутся с зондами, имеющими одинаковый 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 и строятся на графике с двойной логарифмической шкалой.

!Интерактивный Knee Plot (Barcode Rank Plot) для scRNA-Seq: логарифмический график ранжирования баркодов по числу UMI с динамическим порогом отсечения, разделяющим реальные клетки от пустых капель.

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

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

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

Берет пул баркодов с экстремально низким числом UMI (например, ), которые гарантированно являются пустыми каплями, и строит вероятностную модель (распределение Дирихле-мультиномиальное) профиля экспрессии фонового «супа».

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

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

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

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

Общее число молекул / UMI (nCount_RNA). Метрика, тесно коррелирующая с числом генов. Слишком высокие значения также сигнализируют о потенциальных дуплетах.

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

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

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

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

В пространстве транскриптома млекопитающего около 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). Затем для каждого гена вычисляется стандартизированная дисперсия — метрика, показывающая, насколько реальная вариабельность гена превышает ожидаемую для его базового уровня экспрессии.

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

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

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

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

где:

— ковариационная матрица генов, отражающая их совместную изменчивость.

— количество клеток в эксперименте.

— матрица экспрессии, где среднее значение каждого гена сдвинуто к нулю.

— транспонированная матрица экспрессии.

Собственные векторы матрицы задают направления новых осей (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). Практически он означает ожидаемое количество ближайших соседей у каждой клетки, на которые алгоритм обращает внимание.

!Влияние перплексии на t-SNE

Низкая перплексия (например, 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 требует заранее указать количество кластеров (), которое в исследовательском (exploratory) анализе новой ткани неизвестно.

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

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

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

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

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

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

где:

— индекс Жаккара, мера сходства двух множеств (от 0 до 1).

— множество ближайших соседей первой клетки.

— множество ближайших соседей второй клетки.

— количество общих соседей (мощность пересечения множеств).

— общее количество уникальных соседей обеих клеток (мощность объединения множеств).

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

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

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

где:

— модулярность графа (значение от -1 до 1).

— сумма весов всех ребер в графе.

— фактический вес ребра между узлами и .

и — суммы весов всех ребер, присоединенных к узлам и соответственно.

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

— кластеры, к которым отнесены узлы и .

— функция Кронекера: равна 1, если узлы в одном кластере (), и 0 в противном случае.

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

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

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

Главный инструмент биолога при управлении графовой кластеризацией — параметр разрешения (resolution). Математически он внедряется в формулу модулярности как множитель перед штрафом за случайные связи (перед дробью ).

!Влияние параметра resolution на кластеризацию

При низком разрешении (например, 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):

Orchestrating Single-Cell Analysis with Bioconductor (OSCA) — исчерпывающая онлайн-книга (EN) от разработчиков пакетов на R. Содержит глубокий разбор математики отбора признаков (HVG) и графовой кластеризации.

Официальные туториалы пакета Scanpy (Python) на GitHub и ReadTheDocs. В разделе API подробно описана математика под капотом функций sc.tl.pca, sc.tl.umap и sc.tl.leiden.

Официальные виньетки пакета Seurat (R) от Satija Lab. Включают детальные объяснения выбора параметров (таких как resolution и perplexity) на реальных датасетах.

Оригинальные научные публикации (must-read для понимания алгоритмов):

UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction (Leland McInnes et al., 2018) — оригинальная статья, описывающая топологические основы UMAP.

From Louvain to Leiden: guaranteeing well-connected communities (Traag et al., 2019) — статья, математически доказывающая превосходство алгоритма Leiden над Louvain в графовой кластеризации.

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

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

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

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

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

На практике эта стратегия сталкивается с двумя непреодолимыми препятствиями scRNA-Seq: колоссальным размером выборки и экстремальной разреженностью данных (sparsity). В bulk-эксперименте мы сравниваем, например, 3 образца с 3 контролями. В single-cell мы можем сравнивать 5000 клеток одного кластера с 15000 клеток других кластеров. При таких объемах классические параметрические тесты становятся гиперчувствительными: малейшее, биологически незначимое изменение экспрессии гена на доли процента получает астрономически малый -value (например, ). Кроме того, DESeq2 вычислительно тяжел и расчет матрицы на десятки тысяч клеток может занять часы или дни.

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

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

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

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

!Визуализация рангового критерия Уилкоксона

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

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

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

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

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

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

— доля клеток во всех остальных кластерах (фоне), где этот ген детектирован.

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

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

Ген CD3D: , , .

Ген FOXP3: , , .

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

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

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

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

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

!Структура Dot plot для маркерных генов

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

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

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

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

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

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

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

Механика SingleR построена на итеративном вычислении коэффициента ранговой корреляции Спирмена.

На первом шаге алгоритм отбирает гены, которые вариабельны в референсном датасете.

Вычисляется корреляция Спирмена между профилем экспрессии неизвестной клетки и профилем каждого клеточного типа в референсе.

Клетке присваивается предварительная метка того типа, с которым корреляция максимальна.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4. Функциональная аннотация и интерпретация результатов: анализ обогащения 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)

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

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

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

5. Основы 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)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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