Биоинформатика? С R это легко!

Введение 

Эта статья предназначена для тех, кто находится в поиске “точки входа” в область биоинформатики и имеет опыта работы с R (в идеале с использованием tidyverse). Биоинформатика может показаться достаточно пугающей концепцией (в моём случае так и было). Поскольку эта область поистине велика и очень быстро развивается, сходу составить представление о том, что это такое, бывает сложно. Я всегда думал, что биоинформатика — это очень сложная и недоступная для моего понимания сфера, в которой для совершения первых практических шагов сперва потребуется пройти многолетнее нелёгкое обучение. Но, как и в большинстве случаев, оказалось, что для начала требуется не так уж много (при этом для полноценного освоения действительно нужны годы).

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

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

Существует огромное множество разного рода вопросов, равно как и разного рода биологических данных. В этом же примере мы будем выполнять анализ дифференциальной экспрессии генов (DGEA) при помощи данных RNA-Seq (РНК-секвенирования). Коротко говоря, если обратиться к программе колледжа или института, то DGEA отвечает на вопрос, присутствуют ли изменения в уровнях экспрессии генов между разными условиями эксперимента. РНК-секвенирование при этом использует секвенирование следующего поколения для оценки количества РНК в заданном транскриптоне (наборе всех транскрипций РНК). Мы будем получать данные из recount2. По сути, это база данных, администрируемая Университетом Джона Хопкинса, которая взаимодействует с R посредством пакета, а пользователь может запрашивать из неё информацию о подсчётах прочтений в РНК-секвенировании. Подсчёты — это число прочтений в секвенировании, относящихся к конкретной области гена.

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

Настройка

Давайте сначала настроим рабочее пространство, для чего потребуется загрузить в R запрошенные нами из recount библиотеки и данные РНК-секвенирования. Я буду использовать датасет RNA-Seq, относящийся к ишемической болезни сердца, поскольку этот вопрос представляет для меня академический интерес. Кейсы этого датасета относятся к клеткам гладких мышц коронарных артерий с заглушенным геном TCF21. Recount снабжён приятным UI приложения Shiny, через которое можно искать интересующие вас проекты и запрашивать данные, используя ассоциированный ID проекта. Однако в пакете recount также есть функции, позволяющие искать проекты, не обращаясь к интерфейсу приложения. 

library(tidyverse)
library(DESeq2)
library(recount)
library(pheatmap)
library(RColorBrewer)
library(ggrepel)

get_recount_data <- function(projectid) {
  recount::download_study(projectid)
  load(file.path(projectid, "rse_gene.Rdata"),
                envir = .GlobalEnv)
}
get_recount_data("SRP018779")

#Извлекает статус case/control (кейс/контрольные данные)
colData(rse_gene)$case <- colData(rse_gene)$characteristics %>% 
  as.list() %>% 
  map_chr(function(x) ifelse(any(grepl("control", x)),
                             "control",
                             "case"))
#control становится базовым уровнем
colData(rse_gene)$case <- factor(colData(rse_gene)$case, 
                                 levels = c("control", "case"))

В большинстве проектов по анализу данных присутствует 3 основных шага: пробный анализ данных, собственно анализ и визуализация. Прежде чем начать, важно отметить, что для анализа я выбрал DESeq2, хотя существуют и другие пакеты, выполняющие DGEA разными методами. Сначала нам нужно создать объект DESeqDataSet, который мы и будем использовать в анализе. Этот объект, по существу, состоит из трёх компонентов: colData содержит информацию об образцах, rowRanges содержит информацию о genomic ranges (диапазонах геномов) и assay, который содержит матрицу отсчётов. 

Для создания описанного объекта DESeqDataSet мы используем dds <- DESeq2::DESeqDataSet(rse_gene, design = ~ case), где case определяет, является ли образец контрольным или экспериментальным. Теперь давайте перейдём к выполнению первого шага. 

Пробный анализ данных

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

vst_data <- vst(dds, blind = F) #Применение преобразования, стабилизирующего дисперсию

sample_dists <- dist(t(assay(vst_data))) #Вычисление Евклидовых расстояний
sample_dists
sample_dists_mat <- as.matrix(sample_dists)
rownames(sample_dists_mat) <- paste(vst_data$sample, vst_data$case, sep = " - ")
colnames(sample_dists_mat) <- NULL
colors <- colorRampPalette(rev(brewer.pal(9, "Blues")))(255)
pheatmap(sample_dists_mat,
         clustering_distance_rows = sample_dists,
         clustering_distance_cols = sample_dists,
         col = colors, 
         display_numbers = T)
Степень глубины синего определяет степень схожести пар образцов

Ещё один способ оценить схожесть образцов — это использовать график метода главных компонент (PCA), представленный ниже. Не углубляясь в детали, можно сказать, что на этом графике ось x проясняет большую часть отличий между всеми образцами, а ось y вторые по величине различия между ними. Обратите внимание, что эти процентные показатели не суммируются в 100%, потому что есть и другие измерения, которые содержат оставшуюся дисперсию (хотя каждое из них даёт меньше пояснений, чем те, что представлены на графике).

pca_dat <- DESeq2::plotPCA(vst_data, intgroup = "case", returnData = T)
percentvar <- round(100 * attr(pca_dat, "percentVar"))
ggplot(pca_dat, aes(x = PC1, y = PC2, color = case)) +
  geom_point(size = 3) +
  xlab(paste0("PC1: ", percentvar[1], "% variance")) +
  ylab(paste0("PC2: ", percentvar[2], "% variance")) +
  ggtitle("PCA with VST data")
Обратите внимание, как PC1 успешно отделяет кейсы от контрольных образцов. Это говорит о том, что конкретные гены дифференциально экспрессируются. 

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

Теперь, когда мы визуализировали данные, пора перейти к анализу. Используемый нами пакет существенно упрощает задачу, и нам нужно лишь выполнить команду dds ,- DESeq(dds). Вы можете ознакомиться и с другими возможностями этого пакета в R, используя команду помощи: ?DESe. Однако, говоря коротко, выполняет он всего три действия: оценивает факторы размера (контролируя различия в глубине секвенирования, но на данный момент вам об этом задумываться не стоит), оценивает значения дисперсии (тоже не будем углубляться в технические тонкости — можете просто рассматривать дисперсию как меру “распространения”) для каждого гена и выполняет подгонку обобщённой линейной модели. Не пугайтесь, если эти слова показались вам совершенно непонятными. Я придерживаюсь принципа “сначала делать, а потом задавать вопросы”, который в некоторых случаях может быть спорным, хотя на мой взгляд зачастую студенты оказываются настолько перегружены всеми этими теоретическими деталями и техническими аспектами (которые безусловно необходимы для достижения профессионального уровня навыков), что так никогда и не доходят до практической части всего процесса. 

Для каждого гена в датафрейме будет возвращено две основных оценки: baseMean —  среднее из нормализованных значений отсчётов, разделённое на факторы размера, и log2FoldChange — оценка размера эффекта. Например, log2-кратное изменение для 2 означает, что экспрессия гена увеличивается на мультипликативный коэффициент 2²=4. DESeq2 выполняет для каждого гена проверку гипотезы, чтобы понять, достаточно ли доказательств для отклонения нулевой гипотезы, утверждающей, что лечение не оказало на ген влияния, и что причина наблюдаемой разницы между данными, полученными от лечения, и контрольными данными вызвана простой экспериментальной изменчивостью. 

de_data <- DESeq(dds)
res <- results(de_data)
summary(res)

# Из 44567 с ненулевым общим отсчётом прочтений
# Настроенное p-значение < 0.1
# LFC > 0 (повышение)       : 1106, 2.5%
# LFC < 0 (понижение)     : 1251, 2.8%
# Артефакты [1]       : 7210, 16%
# Низкое число прочтений [2]     : 3225, 7.2%
# (среднее число прочтений < 32)

Обратите внимание, что в связи с заглушением TCF21 на уровне FDR в 10% существует множество генов с дифференциальной экспрессией. 

Визуализация

Весьма полезно будет визуализировать распространение оценочных коэффициентов в модели по всем генам. Для этого мы можем использовать MA-график. Log2-кратное изменение (LFC) для конкретного сравнения отражается на оси y, а среднее значение отсчётов, нормализованное по фактору размера, показано на оси x. Каждый ген представлен точкой. Прежде чем создавать MA-график, мы используем функцию lfcShrink, чтобы уменьшить log2-кратные изменения для сравнения наблюдаемых кейсов с контрольными данными. В DESeq2 есть три вида оценки подобного уменьшения. Мы для уменьшения коэффициентов будем использовать метод apeglm, который хорошо справляется с уменьшением зашумленных оценок LFC, давая при этом оценки LFC с низким смещением в случае действительно больших различий. 

dds_shrink <- DESeq2::lfcShrink(de_data, coef = "case_case_vs_control")
res_shrink <- dds_shrink %>%
  as.data.frame() %>%
  rownames_to_column("GeneID") %>%
  left_join(mcols(dds) %>% 
              as_tibble() %>%
              dplyr::select(gene_id, symbol), by = c("GeneID" = "gene_id")) %>%
  dplyr::rename(logFC = log2FoldChange,
                FDR = padj) %>%
  drop_na(pvalue,FDR) %>%
  mutate(GeneID = gsub("\\..*", "", GeneID))
  
  plot_MA <- function(data, names_to_display) {
  cutoff <- sort(data$pvalue)[names_to_display]
  data <- data %>%
    mutate(topgenelabel = ifelse(pvalue <= cutoff, symbol, ""),
           GeneID = gsub("\\..*", "", GeneID))
  ggplot(data, aes(x = log2(baseMean), y = logFC)) + 
    geom_point(aes(col = FDR < 0.05), shape = 20, size = 0.5) + 
    ggrepel::geom_text_repel(aes(label = topgenelabel)) + 
    labs(x = "Mean of Normalized Counts",
         y = "Log Fold Change")
}

plot_MA(res_shrink, 10)
Голубые точки над y=0 показывают гены, существенно повысившие экспрессию, а расположенные ниже y=0 — существенно понизившие. 10 наиболее дифференциально экспрессирующихся генов помечены.

Также всегда полезно выполнять диагностику гистограмм p-значений. 

res_shrink %>%
  filter(baseMean > 1) %>%
  ggplot(aes(x = pvalue)) + 
  geom_histogram(col = "white")
Это нормальная гистограмма с пиком возле 0, показывающим насыщенность дифференциально экспрессирующимися генами. Ничто не указывает на наличие отклонений

Заключение

Итак, чего мы достигли? Мы произвели наш первый анализ дифференциальной экспрессии генов и выявили гены-кандидаты, которые активируются/подавляются, когда заглушен ген TCF21 в клетках мышцы коронарной артерии. Неплохо для начала!

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

Читайте также:

Читайте нас в Telegram, VK и Яндекс.Дзен


Перевод статьи Lathan Liou: How to Start Learning Bioinformatics and Not Get Intimidated (With R)