Методы снижения сложности с данными об экспрессии генов

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

Об остром миелоидном лейкозе: здесь

Об анализе медицинских изображений ИИ: здесь

Об искусственном интеллекте на биологических данных: здесь

Набор данных доступен здесь: набор данных

Набор данных, используемый в этом руководстве, получен от Warnat-Herresthal et al. (1) Они собрали и повторно проанализировали множество наборов данных по лейкемии. В своей статье они разделились на три набора данных (в соответствии с платформой микрочипов). Они включили в базу данных только человеческие образцы трех конкретных платформ микрочипов. В этом уроке для удобства я буду использовать первый набор данных.

Техника микрочипов была разработана в начале 90-х с целью ускорения разработки лекарств. По сути, микрочип — это метод фотолитографии, при котором на чип помещаются тысячи ДНК-зондов. Комплементарная цепь ДНК связывается с ДНК-зондом, и это связывание измеряется. Учитывая возможность одновременного измерения тысяч генов, было потрачено много усилий на разработку биоинформационных методов для анализа данных.

Микрочип позволяет ученым измерять экспрессию мРНК тысяч генов за один раз. В большинстве случаев набор данных состоит из небольшого количества наблюдений и большого количества признаков. Например, набор данных показывает около 2000 образцов (наблюдения, m) и около 14000 генов (признаки, n). В среднем, хотя наборы данных составляют около сотни образцов, они содержат тысячи генов, причем обычно m намного превышает n. Это привело к необходимости методов уменьшения размеров и других инструментов для извлечения важных знаний, скрытых в данных (2). Действительно, этот дисбаланс (m x n) влияет на модель машинного обучения, обученную на наборе данных микрочипа. Многие гены плохо выражены (инвариантны в наборе данных) или избыточны, учитывая, что переменные внутри модели просто вычислительно затратны без какой-либо выгоды. Методы размерного уменьшения могут быть полезны для выбора признаков (найти соответствующие признаки в наборе данных, отбрасывая другие), кластеризации и быстрого взгляда на данные.

Я предположил, что вы уже импортировали набор данных на диск Google.

%reload_ext autoreload
%autoreload 2
%matplotlib inline
from google.colab import drive
drive.mount(‘/content/gdrive’, force_remount=True)

импортировать необходимые библиотеки и набор данных:

#import necessary library
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import umap
#dataset
data = pd.read_table("/content/gdrive/My Drive/aml/201028_GSE122505_Leukemia_clean.txt", sep = "\t")

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

#table of the disease
data.disease.value_counts()

В наборе данных преобладают две категории: острый миелоидный лейкоз (ОМЛ) и острый лимфоидный лейкоз (ОЛЛ). Для удобства визуализации мы сгруппировали некоторые категории.

#removing some disease type
data["disease"] = np.where(data["disease"] == "Diabetes_Type_I" , "Diabetes", data["disease"])
data["disease"] = np.where(data["disease"] == "Diabetes_Type_II" , "Diabetes", data["disease"])
other = ['CML','clinically_isolated_syndrome', 'MDS', 'DS_transient_myeloproliferative_disorder']
data = data[~data.disease.isin(other)]
data.shape
target = data["disease"]
df = data.drop("disease", 1)
df = df.drop("GSM", 1)
df = df.drop("FAB", 1)
df.shape
target.value_counts()

Мы также отфильтровываем функции с низкой дисперсией и масштабируем оставшиеся функции.

df = df.drop(df.var()[(df.var() < 0.3)].index, axis=1)
from scipy.stats import zscore
df = df.apply(zscore)
df.shape

Мы также отфильтровываем функции с низкой дисперсией и масштабируем оставшиеся функции.

df = df.drop(df.var()[(df.var() < 0.3)].index, axis=1)
from scipy.stats import zscore
df = df.apply(zscore)
df.shape

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

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

from sklearn.decomposition import PCA
pca = PCA(n_components=2)
X = pca.fit(df).transform(df)
print(pca.explained_variance_ratio_)

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

from sklearn.preprocessing import LabelEncoder
le = LabelEncoder()
le.fit(target)
y_lan = le.transform(target)
pca_df = pd.DataFrame(columns = [“x”, “y”, “name”, “label”])
pca_df[“PCA1”] = X[:, 0]
pca_df[“PCA2”] = X[:, 1]
pca_df[“Disease”] = target
pca_df[“label”] = y_lan
sns.set(style=”whitegrid”, palette=”muted”)
#sns.set_theme(style=”whitegrid”)
ax = sns.scatterplot(x=”PCA1", y=”PCA2", hue=”Disease”, data=pca_df)
# Put the legend out of the figure
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
#saving the graph in png or jpeg
#plt.savefig(“GSE122505_Leukemia_PCA.pdf”, dpi = 300)
#plt.savefig(“GSE122505_Leukemia_PCA.png”)
#pca_df.to_csv(“GSE122505_Leukemia_PCA.csv”)

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

https://distill.pub/2016/misread-tsne/

Как правило, вы можете следовать, чтобы установить гиперпараметры:

• недоумение должно быть около N^(1/2), где N — количество выборок. Недоумение по умолчанию часто считается равным 30, а другие исследователи предлагают использовать 1% от размера выборки как максимальное недоумение. Недоумение — один из основных параметров, его можно рассматривать как баланс между сохранением глобальной и локальной структуры.
• Количество итераций, в среднем 1000–2000 взаимодействий, должно быть достаточным. Небольшое количество взаимодействий не позволит сформировать кластер, а слишком большое количество итераций слишком долго для вычислений. Эмпирический закон из наблюдения за графиком t-SNE предполагает, что когда наибольшее расстояние между точками данных составляет около 100, алгоритм достигает сходимости.
• Скорость обучения находится в диапазоне от 10 до 1000, слишком высокая и данные слишком разбросаны, слишком низкая и данные слишком сжаты. Скорость обучения по умолчанию во многих библиотеках установлена ​​​​на 100 или 200, некоторые авторы предлагают установить скорость обучения как N/12, если результаты превышают 200.

from sklearn.manifold import TSNE
tsne = TSNE(n_components=2, perplexity=48.0, learning_rate=200.0, n_iter=2000 )
X = tsne.fit_transform(df, y_lan)
tsne_df = pd.DataFrame(columns = ["x", "y", "name", "label"])
tsne_df["tSNE1"] = X[:, 0]
tsne_df["tSNE2"] = X[:, 1]
tsne_df["Disease"] = target
tsne_df["label"] = y_lan

sns.set(style="whitegrid", palette="muted")
ax = sns.scatterplot(x="tSNE1", y="tSNE2", hue="Disease",  data=tsne_df)
# Put the legend out of the figure
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
#plt.savefig("GSE122505_Leukemia_tSNE.pdf", dpi = 300)
#plt.savefig("GSE122505_Leukemia_tSNE.png")
#tsne_df.to_csv("GSE122505_Leukemia_tSNE.csv")

t-SNE может иметь трудности с таким высоким числом размерности, поэтому многие исследователи предлагают уменьшить размерность. Например, вы можете получить это с помощью PCA перед применением t-SNE. В этом случае мы используем PCA и сохраняем первые 50 компонентов. На этих 50 компонентах мы запускаем t-SNE.

pca = PCA(n_components=50)
X_pca = pca.fit(df).transform(df)

tsne = TSNE(n_components=2, perplexity=48.0, learning_rate=200.0, n_iter=2000 )
X = tsne.fit_transform(X_pca, y_lan)
tsne_df = pd.DataFrame(columns = ["x", "y", "name", "label"])
tsne_df["tSNE1"] = X[:, 0]
tsne_df["tSNE2"] = X[:, 1]
tsne_df["Disease"] = target
tsne_df["label"] = y_lan

sns.set(style="whitegrid", palette="muted")
ax = sns.scatterplot(x="tSNE1", y="tSNE2", hue="Disease",  data=tsne_df)
# Put the legend out of the figure
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
#plt.savefig("GSE122505_Leukemia_PCA_tSNE.pdf", dpi = 300)
#plt.savefig("GSE122505_Leukemia_PCA_tSNE.png")
#tsne_df.to_csv("GSE122505_Leukemia_PCA_tSNE.csv")

Различия между PCA и t-SNE:

• PCA - это метод линейного измерения, в то время как t-SNE не является линейным.
• PCA пытается сохранить глобальную структуру данных, в то время как t-SNE пытается сохранить локальную структуру.
• PCA не требует выбора гиперпараметров, а t-SNE — да.
• t-SNE может обрабатывать выбросы, в то время как PCA сильно подвержен влиянию.
• PCA — это детерминированный алгоритм, а t-SNE — рандомизированный.
• PCA вращает векторы, сохраняя дисперсию. Вместо этого t-SNE минимизирует расстояние между двумя точками на гауссовой кривой.

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

Визуализация t-SNE действительно в среднем превосходит визуализацию PCA. Более того, он предпочтительнее транскриптомных данных, поскольку он пытается сохранить локальные структуры. t-SNE генерирует график с разными кластерами, которые часто хорошо согласуются с кластерами, полученными другим алгоритмом кластеризации (3). Однако у самого t-SNE есть некоторые ограничения: он не сохраняет глобальные структуры (межкластерные отношения), требует больших вычислительных ресурсов и не очень осмыслен при работе с очень большими наборами данных (4).

В последние годы новый метод широко используется в контексте визуализации омических данных (часто для секвенирования РНК одной клетки). Утверждается, что равномерное приближение и проекция многообразия (UMAP) превосходит t-SNE в сохранении локальной и глобальной структуры данных и в то же время является более быстрым алгоритмом (4,5). Интересно, что и t-SNE, и UMAP демонстрируют сильную локальную кластеризацию, но аналогичные кластеры в UMAP ближе (глобальная структура). Например, в моде MNIST похожие изображения (например, свитера или пальто) группируются с помощью обоих алгоритмов, но в UMAP аналогичные категории размещаются совместно (например, свитера, рубашки). Это красиво показано здесь.

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

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

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

Однако UMAP также является недетерминированным алгоритмом, и интерпретируемость оси и расстояния не имеет значения (в отличие от PCA). Вам может понадобиться больше, чем график, поскольку стохастическое использование одних и тех же значений параметров приведет к разной визуализации. Кроме того, выбор значения параметра влияет на визуализацию, поэтому лучше попробовать другую комбинацию. Кластеризация исходных данных ограничена проклятием размерности (упомянутый выше дисбаланс m x n), в то время как на t-SNE она довольно сложна, но ее можно легко сделать с результатами UMAP и PCA (о чем будет подробно рассказано в следующем уроке).

import umap
reducer = umap.UMAP(n_neighbors =  100, min_dist= 0.2, metric ="euclidean")
X = reducer.fit_transform(df)
umap_df = pd.DataFrame(columns = ["x", "y", "name", "label"])
umap_df["UMAP1"] = X[:, 0]
umap_df["UMAP2"] = X[:, 1]
umap_df["Disease"] = target
umap_df["label"] = y_lan

sns.set(style="whitegrid", palette="muted")
ax = sns.scatterplot(x="UMAP1", y="UMAP2", hue="Disease",  data=umap_df)
# Put the legend out of the figure
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
#plt.savefig("GSE122505_Leukemia_UMAP.pdf", dpi = 300)
#plt.savefig("GSE122505_Leukemia_UMAP.png")
#ctrl_df.to_csv("GSE122505_Leukemia_UMAP.csv")

Точно так же мы можем сначала использовать PCA, а затем запустить UMAP.

pca = PCA(n_components=50)
X_pca = pca.fit(df).transform(df)

reducer = umap.UMAP(n_neighbors =  100, min_dist= 0.2, metric ="euclidean")
X = reducer.fit_transform(X_pca)

umap_df = pd.DataFrame(columns = ["x", "y", "name", "label"])
umap_df["UMAP1"] = X[:, 0]
umap_df["UMAP2"] = X[:, 1]
umap_df["Disease"] = target
umap_df["label"] = y_lan

sns.set(style="whitegrid", palette="muted")
ax = sns.scatterplot(x="UMAP1", y="UMAP2", hue="Disease",  data=umap_df)
# Put the legend out of the figure
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
#plt.savefig("GSE122505_Leukemia_PCA_UMAP.pdf", dpi = 300)
#plt.savefig("GSE122505_Leukemia_PCA_UMAP.png")
#ctrl_df.to_csv("GSE122505_Leukemia_PCA_UMAP.csv")

Результаты схожи в обоих случаях.

Затем мы можем использовать UMAP для различного исследовательского анализа. В проекции UMAP мы можем выделить конкретный ген, в данном случае ANXA1.

umap_df["ANXA1"] = df["ANXA1"]
sns.set(style="whitegrid", palette="muted")
ax = sns.scatterplot(x="UMAP1", y="UMAP2", hue="ANXA1",  data=umap_df)
# Put the legend out of the figure
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)

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

Избранная библиография

1. Musheer R, Verma CK, Srivastava N. Методы уменьшения размерности для данных микрочипов: обзор. АИМС Биоинженерия. 2017;4:179–97.

2. Варнат-Херресталь С., Перракис К., Ташлер Б., Беккер М., Басслер К., Бейер М. и др. Масштабируемое прогнозирование острого миелоидного лейкоза с использованием многомерного машинного обучения и транскриптомики крови. iScience [Интернет]. 2019 [цитировано 26 октября 2020 г.]; 23. Доступно по адресу: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6992905/

3. Кобак Д., Беренс П. Искусство использования t-SNE для транскриптомики одиночных клеток. Связь с природой. Издательская группа «Природа»; 2019;10:5416.

4. Becht E, McInnes L, Healy J, Dutertre C-A, Kwok IWH, Ng LG, et al. Уменьшение размерности для визуализации данных одной ячейки с использованием UMAP. Природная биотехнология. Издательская группа «Природа»; 2019;37:38–44.

5. Макиннес Л., Хили Дж., Мелвилл Дж. UMAP: Единообразная аппроксимация многообразия и проекция для уменьшения размерности. arXiv:180203426 [cs, stat] [Интернет]. 2020 [цитировано 28 октября 2020 г.]; Доступно по адресу: http://arxiv.org/abs/1802.03426

Доступность кода

Код урока и все изображения вы найдете здесь. Вы можете использовать предоставленный блокнот Colab без необходимости установки каких-либо пакетов на свой компьютер.

если вам было интересно:

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