Биоинформатика #20: Предсказание экспрессии генов миелодиспластических синдромов с помощью машинного обучения…

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

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

Анализ предсказания генов

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

Я ищу несколько статей и нахожу интересную статью исследовательской группы Лю еще в 2019 году, когда они написали статью о прогнозировании экспрессии генов с использованием обобщенной линейной модели (GLM). Когда я прочитал их статью, и они были процитированы более сотни раз, я пришел к выводу, что предсказание генов — это хорошее приложение, которое необходимо в биомедицинской науке. Модель прогнозирования можно использовать для проверки взаимосвязи между геном или комбинацией генов (входные данные) и желаемым результатом фенотипа, таким как промоторная область или даже экспрессия другого гена (выходные данные).

Одной из моделей машинного обучения, которую я хочу применить, является линейная модель, которая требует непрерывных входных данных и непрерывных выходных данных. Поскольку исследование Лю было применено в системе линейных моделей, мне интересно резюмировать их исследовательскую модель, но я буду использовать другой набор данных. Набор данных, который я нашел, — это миелодиспластические синдромы, в которых для проверки экспрессии генов использовался микрочип (Kim, et al., 2015).

Миелодиспластические синдромы

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

Я проверяю их данные на NCBI GEO, чтобы получить данные микроматрицы (входные данные), а также проверяю данные о тромбоцитах пациента (PLT) (выходные данные). Моя цель состоит в том, я хочу знать, можно ли использовать 20 самых значимых генов для предсказания тромбоцитов пациентов. Будем надеяться, что если модель машинного обучения (ML) будет хорошей, то 20 наиболее значимых генов в этом наборе данных можно будет использовать в качестве биомаркеров для обнаружения клеток крови тромбоцитов с миелодиспластическим синдромом.

Процедура анализа

Анализ начинается с анализа генов дифференциальной экспрессии, потому что мы хотим найти 20 наиболее значимых генов в данных микрочипа с использованием пациентов с миелодиспластическим синдромом по сравнению с нормальными пациентами. Затем мы должны прочитать исходную статью(Kim, et al., 2015), чтобы найти количество тромбоцитов у каждого пациента в качестве выходных данных. Далее мы объединяем их в одну таблицу и анализируем ее с помощью пакета scikit от python. Для приложения пакета scikit в инструкции python я отдаю должное статье Роберта Таса Джона, в которой объясняется применение линейной модели регрессии ML.

1. Анализ генов дифференциальной экспрессии

Для анализа генов дифференциальной экспрессии (DEG) мы будем использовать GEO2R, который позволит автоматически анализировать данные микрочипов (Barret, et al., 2012). Плюсом GEO2R является то, что они также предоставляют исходный код R, который мы можем использовать для дополнительной модификации (рис. 1). В этом случае я выбираю выборку и классифицирую ее на «нормальных» и «больных». Затем я изменяю код, чтобы получить необработанные данные об экспрессии генов, которые можно использовать для входных данных машинного обучения линейной модели регрессии.

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

# load series and platform data from GEO

gset <- getGEO("GSE61853", GSEMatrix =TRUE, AnnotGPL=TRUE)
if (length(gset) > 1) idx <- grep("GPL10558", attr(gset, "names")) else idx <- 1
gset <- gset[[idx]]

# make proper column names to match toptable 
fvarLabels(gset) <- make.names(fvarLabels(gset))

# group membership for all samples
gsms <- "11111110000000"
sml <- strsplit(gsms, split="")[[1]]

# log2 transformation
ex <- exprs(gset)
qx <- as.numeric(quantile(ex, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
LogC <- (qx[5] > 100) ||
  (qx[6]-qx[1] > 50 && qx[2] > 0)
if (LogC) { ex[which(ex <= 0)] <- NaN
exprs(gset) <- log2(ex) }

# assign samples to groups and set up design matrix
gs <- factor(sml)
groups <- make.names(c("SYNDROMES","CONTROL"))
levels(gs) <- groups
gset$group <- gs
design <- model.matrix(~group + 0, gset)
colnames(design) <- levels(gs)

fit <- lmFit(gset, design)  # fit linear model

# set up contrasts of interest and recalculate model coefficients
cts <- paste(groups[1], groups[2], sep="-")
cont.matrix <- makeContrasts(contrasts=cts, levels=design)
fit2 <- contrasts.fit(fit, cont.matrix)

# compute statistics and table of top significant genes
fit2 <- eBayes(fit2, 0.01)
tT <- topTable(fit2, adjust="fdr", sort.by="B", number=250)

tT <- subset(tT, select=c("ID","adj.P.Val","P.Value","t","B","logFC","Gene.symbol"))
tT[tT == ""] <- NA 
tT <- tT[!is.na(tT$Gene.symbol),]
write.table(tT, file=stdout(), row.names=F, sep="\t")


#Top 20 significant DEGs
top_20_genes <- tT[1:20,]
top_20_genes_code <- rownames(top_20_genes) #extract top 20 Illumina code for the raw gene expression extraction

top_20_genes_expression <- ex[rownames(ex) %in% top_20_genes_code,]
top_20_genes_code %in% rownames(top_20_genes_expression) #Check either the gene expression ID is extracted correctly

2. Извлеките информацию из исходной бумаги.

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

3. Модель линейной регрессии машинного обучения

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

#Linear regression modelling

steps = [
    ('scalar', StandardScaler()),
    ('poly', PolynomialFeatures(degree=2)),
    ('model', LinearRegression())
]

pipeline = Pipeline(steps)

pipeline.fit(X_training, Y_train)

print('Training score: {}'.format(pipeline.score(X_training, Y_train)))
print('Test score: {}'.format(pipeline.score(X_test, Y_test)))

В этом случае мы должны построить конвейер, чтобы упростить управление процессом машинного обучения (John, 2018). В этом случае я выбираю полиномиальные признаки, потому что мы хотим знать вклад каждого гена из 20 генов в модель линейной регрессии. Результат довольно плохой, потому что оценка за обучение равна 100%, а оценка за тест — минус. Это означает, что с входными данными что-то не так, причину я объясню позже. На данный момент мы можем использовать регуляризацию для повышения результатов теста.

4. Регуляризация

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

#Regularization using lasso regression
steps = [
    ('scalar', StandardScaler()),
    ('poly', PolynomialFeatures(degree=2)),
    ('model', Lasso(alpha=0.2, fit_intercept=True))
]

lasso_pipe = Pipeline(steps)

lasso_pipe.fit(X_train, y_train)

print('Training score: {}'.format(lasso_pipe.score(X_training, Y_train)))
print('Test score: {}'.format(lasso_pipe.score(X_test, Y_test)))

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