За последнее десятилетие количество данных в биологии резко возросло, и, похоже, для этого нет предела, поскольку у нас есть все больше и больше секвенированных геномов, более глубокая экспрессия РНК и наборы данных эпигенетики. Отдельно отметим, что у нас есть блестящие примеры широкого использования машинного обучения для предсказания структуры белков и открытия лекарств. Теперь мы также можем изучать трехмерную организацию генома, которая позже может быть использована для понимания геномных перестроек и последующих заболеваний.
Но есть одна относительно новая область биоинформатики, в которой классические инструменты машинного обучения особенно полезны - секвенирование одноклеточной РНК.
Прежде чем перейти к данным по отдельным клеткам, давайте кратко рассмотрим, что такое секвенирование РНК и что такое РНК в целом.
Некоторые из вас, возможно, слышали о центральной догме молекулярной биологии, например, о преобразовании информации из ДНК в РНК и из РНК в белки. Белки - это молекулярные машины, отвечающие за все в клетке. К сожалению, очень сложно количественно определить точное количество каждого белка в клетке, но мы легко можем секвенировать их соответствующую РНК.
Как указано в названии, секвенирование РНК - это группа методов для количественного определения количества РНК в образце. Причина, по которой мы хотим дать количественную оценку, заключается в том, что разные клетки могут экспрессировать разные гены на разных уровнях. Из-за этой разницы, клетки обладают чрезвычайно разными свойствами и функциями, и нам интересно, что их вызвало. Типичными примерами использования секвенирования РНК является сравнение «активности» гена (количества РНК) в двух группах, например, лечение и контроль, где у нас есть N контрольных образцов, N образцов лечения, и мы используем статистические методы для оценки того, какие гены почтительно экспрессируются (= имеют разные средние значения в популяциях) между нашими группами (это можно интерпретировать напрямую как тест A \ B для тех, кто знаком с ними). Это очень полезно, когда мы пытаемся понять, оказывает ли, например, наше лечение какой-либо эффект. Другой пример - наша попытка понять влияние различных мутаций геномной последовательности на уровень экспрессии генов.
Однако у этих методов есть свои ограничения. Поскольку мы извлекаем РНК из группы клеток, мы теряем всю информацию о различиях между ними и можем получить только среднее значение экспрессии каждого гена. В некоторых случаях это так же полезно, как знать среднюю температуру всех пациентов в больнице. Это может быть число, но имеет ли он смысл? Это особенно важно для исследований рака, поскольку было показано, что опухоль неоднородна, но на самом деле может состоять из клетки с разными морфологическими, фенотипическими профилями и профилями экспрессии, поэтому тесты на лекарственную чувствительность могут быть не на 100% точными. Однако это не означает, что массовое секвенирование РНК бесполезно, нам просто нужно понимать его ограничения.
так что нам делать? Что ж, не углубляясь в детали молекулярной биологии, мы можем попытаться сохранить информацию о каждой клетке, поэтому каждый раз, когда мы секвенируем молекулу РНК, мы также секвенируем небольшую техническую последовательность, которая уникальна для каждой клетки. Если мы дополнительно добавим небольшую техническую последовательность к каждой уникальной молекуле РНК, мы также сможем точно определить реальное количество РНК в каждой клетке. Таким образом, для каждой молекулы мы секвенируем последовательность, мы знаем идентификатор ячейки (штрих-код), из которой она произошла, и ее собственный уникальный идентификатор (в биоинформатике мы называем это уникальным молекулярным идентификатором, UMI), чтобы мы могли количественно определить точное количество РНК в каждой ячейке. , поэтому сохраняю всю возможную информацию! В конце концов, в конце эксперимента у нас будет таблица P по генам, где P [i, j] - это экспрессия гена i в ячейке j.
Но, как и у любой техники, есть проблемы и ограничения. И одна из таких проблем, самая важная для нас, - это выпадение учащихся. Из-за стохастической природы экспрессии генов и технических ограничений от каждой клетки мы получаем только подмножество генов, подмножество данных, которые мы теоретически можем получить из них. Мы наблюдаем много нулевой экспрессии различных генов в наших клетках, и, к сожалению, мы не можем сказать, является ли один конкретный «ноль» потому, что этот ген просто не экспрессируется в клетке или он просто не попал в нашу подгруппу - мы называем такие события «выпадением». (Некоторые люди могут заметить, что это очень похоже на системы рекомендаций, где у нас также довольно скудные данные, поскольку один пользователь не может оценить все возможные фильмы \ продукты)
Теперь, когда мы рассмотрим все детали, необходимые для понимания наших данных, как с биологической, так и с технической точки зрения, мы можем резюмировать их с точки зрения задачи машинного обучения.
Итак, наши данные обладают следующими свойствами:
Типичное количество признаков (генов) составляет ~ 20 000–30 000. (Здесь, однако, мне не нужно обращать внимание на то, что мы обычно удаляем гены с низкой экспрессией и работаем примерно с ~ 5 000–7 000 генов… что по-прежнему много). Количество образцов может варьироваться от 5000 ячеек (обычное количество от одного прохода конвейера, который мы используем в лаборатории) до 400000 ячеек и даже 1300000 ячеек. Как упоминалось выше, из-за отсева и характера выборки данные скудны. Еще один важный момент для алгоритмов кластеризации - это разные масштабы наших характеристик (ген A может быть экспрессирован в 100 раз выше, чем ген B, но это не значит, что ген A важнее, чем ген B) и масштабы выборки. У нас может быть разное количество UMI (сумма экспрессии) на клетку, но это может быть биологическая вариация (некоторые клетки в течение своей жизни просто имеют меньше мРНК) или техническая (две клетки могут иметь точно такое же количество мРНК, но из one мы получаем больше данных из-за технических вариаций). Поэтому правильная нормализация имеет решающее значение.
Типичный вопрос, который задают биологи, обычно задают всем K-клеткам, сколько состояний клеток (кластеров клеток) у нас есть и каковы их гены-маркеры, которые мы определяем как статистически значимые характеристики определенного кластера. Учитывая, что у нас редко есть какие-либо ярлыки из ячеек, мы можем сказать, что это типичный вопрос для обучения без учителя, когда у нас нет предварительных знаний о количестве кластеров, и кластеры могут быть разного размера , форма, плотность и могут быть связаны друг с другом (например, клеточный цикл добавляет дополнительную вариацию между близкими состояниями клеток). Принимая во внимание компонент выпадения, который добавляет дополнительный шум к данным и может увеличивать расстояние между близкими ячейками, эта задача становится все более сложной.
Другой важной проблемой здесь является высокая размерность наших данных из-за так называемого «проклятия размерности». Поскольку мы работаем как минимум с 5000 генами, потребуется 2⁵⁰⁰⁰ клеток, чтобы покрыть все точки в этом пространстве, даже если бинаризовать экспрессию генов. Очевидно, что мы не можем упорядочить такое количество ячеек, поэтому нам нужно сначала уменьшить количество измерений.
Типичный конвейер для этого - поиск сильно изменчивых функций (например, с использованием коэффициента вариации), а затем использование PCA для уменьшения количества компонентов до некоторого разумного количества и, кроме того, для визуализации, использование tSNE или UMAP или даже более сложные методы, такие как вариационные автоэнкодеры и даже GAN (однако это может работать только с большими наборами данных, и такие модели сложно перенести из одного набора данных в другой).
Следующий шаг в анализе - кластеризация. В настоящее время широко используются по крайней мере 15 различных инструментов для кластеризации последовательностей одиночных клеток РНК с использованием различных механизмов кластеризации, от простых K-средств с несколькими инициализациями до dbscan и агломеративной кластеризации. Однако, поскольку у нас обычно отсутствуют какие-либо метки для ячеек, мы можем попытаться оценить надежность наших кластеров на наборах данных с известными метками или с помощью различных метрик для оценки производительности кластера, таких как BIC, AIC или метод силуэта, или мы можем сделать умнее - попросите помощи у биологов! Мы можем использовать знания в предметной области, поскольку для многих важных генов мы знаем характер их экспрессии и, используя это, мы можем попытаться проанализировать и идентифицировать некоторые известные типы клеток. Это не всегда возможно, но это тема для дальнейших статей.
Давайте построим простой конвейер для анализа РНК одиночных клеток. Мы проанализируем общедоступный набор данных мононуклеарных клеток периферической крови (PBMC) от 10X Genomics, который состоит из ~ 8000 клеток, который доступен здесь
Далее мы загружаем наши данные. Матрица довольно большая с большим количеством нулей, поэтому хранится в разреженном формате.
Давайте сделаем EDA, чтобы понять, с чем мы работаем
Мы видим, что у нас много генов с чрезвычайно низкой экспрессией, а некоторые (возможно, даже один с чрезвычайно высокой экспрессией) - они могут быть выбросами. С другой стороны, мы также знаем, что некоторые факторы транскрипции (гены, регулирующие экспрессию других генов) могут экспрессироваться на очень низком уровне, и фильтрация наших генов может привести к потере некоторой структуры.
В то же время работать с 33555 функциями довольно сложно, поэтому давайте удалим некоторые из них. Для этого мы можем сделать простую вещь. Сначала мы удаляем гены с низкой экспрессией (скажем, 100), строим график CV (коэффициент вариации) и выбираем разумный порог, выбираем гены, которые мы будем использовать для кластеризации.
Давайте определим высоко вариабельные и полезные гены как гены с CV ≥ 10.
Нам снова нужно проверить наши клетки. Может случиться так, что после фильтрации генов мы получим клетки с нулевой экспрессией каких-либо генов. Лучше удалить их и после лог нормализовать наши данные
Следующим нашим шагом будет использование UMAP. Мы будем использовать его, чтобы уменьшить количество измерений до 2, и мы запустим кластеризацию для этого. Это не лучшая идея, поскольку проекция только по двум измерениям может быть не оптимальной и, вероятно, не будет отражать всю сложность набора данных, но у нее есть один плюс - легко визуализировать и понимать результаты различных методов кластеризации. Еще мы можем попробовать оптимизировать здесь - это n_neighbors, min_dist, скорость обучения и другие параметры UMAP.
Обычно кластеризация выполняется после PCA, и мы используем tSNE для проецирования результатов на двумерную плоскость. Однако UMAP лучше, если мы захотим добавить дополнительные точки данных позже. Ну, UMAP в целом лучше.
Мы уже видим здесь некоторую структуру, без оптимизации параметров UMAP и pca и с довольно наивным способом выбора функций для кластеризации. Вопрос только в том, отражает ли эта структура биологию. Но в качестве упражнения мы можем попробовать запустить несколько методов кластеризации.
Просто визуально DBSCAN и SpectralClustering выглядят хорошо, но, как я уже сказал, есть некоторые показатели, которые нужно отслеживать, чтобы найти «оптимальную» кластеризацию. Поэтому здесь мы также можем настроить наши параметры. В конце концов, наш наивный подход можно улучшить за три шага.
- Точный выбор функций
- Настройка параметров для UMAP + pca
- Настройка параметров для кластеризации
Однако главный вопрос здесь, как я уже сказал, отражает ли эта структура истинную биологическую сложность. Чтобы ответить на этот вопрос, нам нужно будет взглянуть на гены-маркеры каждого идентифицированного кластера, чтобы попытаться сопоставить их с известными типами клеток и понять лежащую в основе биологию.
В следующем посте мы сравним наши результаты с часто используемыми конвейерами и попытаемся найти «маркеры» для каждого кластера и идентифицировать их.
