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

РНК-полимераза считывает ДНК и синтезирует РНК в ходе транскрипции.

Что дает секвенирование РНК?

На Habr уже присутствуют статьи, которые рассказывают о том, какие инструменты применяются в биоинформатике для решения научных задач в биологии. Так, например, есть достаточно интересная серия статей «Практическая биоинформатика». Здесь хотелось бы затронуть целенаправленно такую область, как транскриптомика на языке R. Мы попробуем максимально наглядно и с пояснениями разобрать что такое это за тип анализа, как получают данные и для каких целей он применяется.

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

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

Транскриптомика делится на два больших раздела — bulk RNA-seq и single-cell RNA-seq. В названии отражена и суть данных методов: bulk RNA-seq - это секвенирование ткани целиком, тогда как single-cell RNA-seq - каждой отдельной клетки и сохранение информации об экспрессии. Здесь будет подробно освещен пайплайн работы с bulk RNA-seq.

Что это вообще такое? Фактически, мы берем для секвенирования весь образец ткани, всю совокупность клеток, молекулы РНК записываются в строковом виде фрагментов по 100-200 нуклеотидов. Затем мы собираем из них целостные РНК и считаем экземпляры повторяющихся строк. Потом объединяем все образцы с посчитанными генами в одну матрицу, количество столбцов в которой равно количеству образцов, а строк — генам. У человека будет около 40 тыс. строк в матрице.

В рамках транскриптомики измеряется дифференциальная экспрессия генов. То есть находится разница в количестве молекул РНК у контроля и эксперимента, или у одного образца по сравнению с другим, если, например, анализируется динамика заболевания. При этом изменение экспрессии описывается логарифмической функцией, из чего следует, что, при его отрицательном значении, образец показывает снижение синтеза соответствующего гена, а при положительном — повышение. При этом bulk RNA-seq дает информацию лишь о ткани или образце в среднем, хотя такой анализ можно провести даже на слабом ноутбуке, тогда как single-cell охватывает данные о каждой индивидуальной клетке, которых могут быть тысячи. Количество информации б полученной в ходе single-cell анализа также в тысячи раз больше, и потому такой анализ предпочтительно проводить на сервере, либо, если данных не очень много (например, несколько сотен клеток) - на мощном компьютере.

Кратко резюмируем необходимую терминологию

Хоть мы и будем говорить в основном о том, «как делать?» с точки зрения написания кода и визуализаций, нужно всё же пояснить некоторые биологические термины, чтобы не теряться в ходе повествования. Их немного:

Секвенирование — фактически, чтение геномных данных образца, превращение биологических молекул в строковую информацию;
RNA-seq — как ни странно, секвенирование РНК;
bulkRNA-seq — секвенирование образца ткани целиком, без разделения на отдельные клетки. Менее вычислительно затратен, т.к. у нас есть один экземпляр матрицы экспрессии. Фактически, «среднее по больнице» количество генов для всего образца;
Экспрессия — процесс превращения информации в белки внутри клетки. В нашем случае, речь идет именно о транскрипции, а количественный анализ проводится для РНК. Любой стимул, любое воздействие или процесс регулируется белками. Экспрессия генов, кодирующих их структуру, меняется в ответ на такие воздействия;
Single-cell RNA-seq — то же самое, только для каждой клетки есть информация об экспрессии. Фактически, объем данных bulkRNA-seq умножается на количество клеток в образце;
Дифференциальная экспрессия — логарифмический показатель изменения экспрессии в образце по сравнению с контролем. Может быть как отрицательным (экспрессия снижается), так и положительным;
Прочтения, они же "риды" (от "read") - короткие фрагменты РНК, которые записываются в FASTQ файлах в виде 4-строчных последовательностей символов (в статье далее пример);
Референс - эталон, заранее известный геном исследуемого организма, с которым происходит сравнение, на основе которого из ридов собирается уже секвенированный геном путем сопоставления.

Что используют для анализа экспрессии?

Основным языком в деле анализа данных долгое время был R. Пакеты, которые понадобятся:

Визуализация: gplots, ggplot2, plotly, ComplexHeatmap, EnhancedVolcano, RColorBrewer, circlize;

Работа с информацией непосредственно: org.Hs.eg.db, DESeq2, edgeR, ReactomePA, clusterProfiler.

Для bulkRNA-seq данных обычно требуются совсем небольшие ресурсы — до 16 ГБ ОЗУ хватит для большинства операций. Для scRNA-seq используют сервер или кластер, и пакеты Seurat, Matrix, celldex, SingleR и другие. Если одна матрица bulkRNA-seq - это около 40 тыс. строк (Картинка 2) с десятками столбцов, то даже просто входные данные scRNA-seq - те же данные, умноженные на количество клеток, которых может быть более 100 тыс.

Картинка 2. Пример матрицы экспрессии, фрагмент.

Столбцы — названия образцов, строки - entrez id генов, в ячейках - количество молекул РНК каждого гена для каждого образца.

Объем входной матрицы bulkRNA-seq обычно < 1 МБ, а для scRNA-seq может быть >1 ГБ.

Сравнение происходит между столбцами матрицы. Как правило, несколько из них — контроль, полученный у здоровых людей. Остальные - экспериментальные образцы. Например, ткани, подверженной болезни, или какому-либо другому воздействию.

Проводится сравнение с использованием теста Вальда с поправкой на множественные сравнения Бенджамини-Хохберга. Впрочем, это регулируется, вариант поправки может быть задан в аргументе функции DESeq:

DESeq(object, test = c(«Wald», «LRT»), 
      fitType = c(«parametric», «local», «mean», «glmGamPoi»), 
      sfType = c(«ratio», «poscounts», «iterate»), 
      betaPrior, 
      full = design(object), 
      reduced, 
      quiet = FALSE, 
      minReplicatesForReplace = 7, 
      modelMatrixType, 
      useT = FALSE,
      minmu = if (fitType == «glmGamPoi») 1e-06 else 0.5, 
      parallel = FALSE, 
      BPPARAM = bpparam())

Как получают данные? И очень коротко про предобработку…

Для описанного выше необходимы, как ни странно, транскриптомные данные. Данные РНК-секвенирования интересующей ткани. Изначально это миллионы фрагментов молекул РНК в формате FASTQ.

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

@Some RNA-seq sample 1
CGUAGCUGCUCGAUCGAUCGAUCGAUCAGU
+
()(&&*III#II@@+++//%$$#(&A

Таблица кодирования quality score в 4 строчке FASTQ-рида

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

Однако, первым шагом будет очистить данные от низкокачественных прочтений. Секвенатор заранее записывает, какие нуклеотиды удалось считать хорошо, а какие — не очень. Такие инструменты, как FastQC и MultiQC позволяют оценить качество. Кроме того, нужно удалить вспомогательные цепочки - адаптеры. Они стандартные, поэтому можно напрямую использовать заранее известную строковую последовательность из данных при помощи, например, cutadapt. Отфильтровав низкокачественные прочтения, а также отрезав адаптеры, но оставив дупликаты (далее станет ясно, почему), можно переходить к выравниванию - процессу собирания из разрозненных прочтений целостных генов, которые соответствуют информации в базах. Для этого используют референсные последовательности из NCBI, например, генома человека, и инструменты для выравнивания. Это может быть как Kallisto, основанное на ML, так и, например, minimap2, bowtie2 и другие инструменты. Как выглядит выравнивание, можно увидеть на Картинке 3.

Картинка 3. Схема процесса выравнивания рида на референс.

Вкратце, наши риды сопоставляются с референсом. Например, у нас есть геном человека — длинная цепь из нуклеотидов. Риды накладываются на него, благодаря чему можно понять, в каком порядке они изначально были расположены в геноме, и собрать такие же целостные гены, как в геноме, но уже с нашего образца.

На выходе получаются файлы формата BAM, в которых те же последовательности, но теперь заметно более длинные. Это и есть гены. Чтобы собрать их в матрицу, их нужно просто сосчитать. Так как мы считаем количество экземпляров повторяющихся строк, нельзя удалять дубликаты, как это обычно делается — нам нужно количественно оценить экспрессию, а значить, по определению, каких-то генов было экспрессировано больше, а каких-то - меньше. Подсчет можно производить, например, при помощи Kallisto quant или Salmon quant. В итоге, для каждого образца получается вектор, а точнее, матрица с одним столбцом. Названием столбца является название образца, строки - ген. Понятно, что так как секвенируют в одном эксперименте обычно одни и те же организмы, количество генов у них одинаковое, а значит, для каждого образца длина вектора названий генов будет одинаковой. Их совмещают, таким образом, получается матрица, где столбцами являются образцы, строками - гены. (Картинка 2) А в ячейках представлены количества РНК-молекул каждого гена для каждого образца. Так и получаются матрицы экспрессии. Более подробно о предобработке - в этой статье на Habr. Подробный же туториал обработки данных секвенирования РНК можно найти здесь.

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

Итак, закончили с пояснениями процесса, перейдем к коду. Для нашего примера, возьмем данные секвенирования у пациентов с нарушением обонятельной функции после коронавируса. Сравнивать будем образцы у переболевших (это «эксперимент», COVID-19 Positive) со здоровыми (»контроль», COVID-19 Negative)

# Считывание матрицы экспрессии
df <- read.csv('raw_counts.tsv', sep = '\t')
names <- colnames(df[,-1])

# Метаданные эксперимента
meta <- read.csv('metadata.txt', sep = ',')

# Переименование столбцов для удобства
colnames(meta) <- c('sample', 'status')
metastatus)
genes <- df$GeneID
rownames(df) <- genes
df$GeneID <- NULL
genes <- as.character(genes)

# Создание объекта для дальнейшего расчета дифференциальной экспрессии
dds <- DESeqDataSetFromMatrix(countData = df,
                              colData = meta,
                              design = ~ status)

# Используем общепринятые наименования генов
symbols <- mapIds(org.Hs.eg.db, keys = genes,
                  column = c('SYMBOL'), 
                  keytype = 'ENTREZID')

# Фильтруем пропущенные значения
dds <- dds[!is.na(rownames(dds)),]

# Запуск анализа
dds <- DESeq(dds)
res <- results(dds)

# Сохранение значимых дифференциально-экспрессирующихся генов для визуализации
significant_genes <- res[!is.na(respadj < 0.05 & 
                           !is.na(res$log2FoldChange) & 
                           abs(res$log2FoldChange) > 1, ]

Части со считыванием данных и сохранением в переменную вряд ли нуждаются в пояснении. То, что идет после считывания — создания объекта «dds» для расчета дифференциальной экспрессии, как это предусмотрено в DESeq2. Структура «dds» выглядит так:

Formal class 'DESeqDataSet' [package "DESeq2"] with 8 slots
  ..@ design            :Class 'formula'  language ~status
  .. .. ..- attr(*, ".Environment")= 
  ..@ dispersionFunction:function ()  
  ..@ rowRanges         :Formal class 'CompressedGRangesList' [package "GenomicRanges"] with 5 slots
  .. .. ..@ unlistData     :Formal class 'GRanges' [package "GenomicRanges"] with 7 slots
  .. .. .. .. ..@ seqnames       :Formal class 'Rle' [package "S4Vectors"] with 4 slots
  .. .. .. .. .. .. ..@ values         : Factor w/ 0 levels: 
  .. .. .. .. .. .. ..@ lengths        : int(0) 
  .. .. .. .. .. .. ..@ elementMetadata: NULL
  .. .. .. .. .. .. ..@ metadata       : list()
  .. .. .. .. ..@ ranges         :Formal class 'IRanges' [package "IRanges"] with 6 slots
  .. .. .. .. .. .. ..@ start          : int(0) 
  .. .. .. .. .. .. ..@ width          : int(0) 
  .. .. .. .. .. .. ..@ NAMES          : NULL
  .. .. .. .. .. .. ..@ elementType    : chr "ANY"
  .. .. .. .. .. .. ..@ elementMetadata: NULL
  .. .. .. .. .. .. ..@ metadata       : list()
  .. .. .. .. ..@ strand         :Formal class 'Rle' [package "S4Vectors"] with 4 slots
  .. .. .. .. .. .. ..@ values         : Factor w/ 3 levels "+","-","*": 
  .. .. .. .. .. .. ..@ lengths        : int(0) 
  .. .. .. .. .. .. ..@ elementMetadata: NULL
  .. .. .. .. .. .. ..@ metadata       : list()
  .. .. .. .. ..@ seqinfo        :Formal class 'Seqinfo' [package "GenomeInfoDb"] with 4 slots
  .. .. .. .. .. .. ..@ seqnames   : chr(0) 
  .. .. .. .. .. .. ..@ seqlengths : int(0) 
  .. .. .. .. .. .. ..@ is_circular: logi(0) 
  .. .. .. .. .. .. ..@ genome     : chr(0) 
  .. .. .. .. ..@ elementMetadata:Formal class 'DFrame' [package "S4Vectors"] with 6 slots
  .. .. .. .. .. .. ..@ rownames       : NULL
  .. .. .. .. .. .. ..@ nrows          : int 0
  .. .. .. .. .. .. ..@ elementType    : chr "ANY"
  .. .. .. .. .. .. ..@ elementMetadata: NULL
  .. .. .. .. .. .. ..@ metadata       : list()
  .. .. .. .. .. .. ..@ listData       : Named list()
  .. .. .. .. ..@ elementType    : chr "ANY"
  .. .. .. .. ..@ metadata       : list()
  .. .. ..@ elementMetadata:Formal class 'DFrame' [package "S4Vectors"] with 6 slots
  .. .. .. .. ..@ rownames       : NULL
  .. .. .. .. ..@ nrows          : int 39376
  .. .. .. .. ..@ elementType    : chr "ANY"
  .. .. .. .. ..@ elementMetadata:Formal class 'DFrame' [package "S4Vectors"] with 6 slots
  .. .. .. .. .. .. ..@ rownames       : NULL
  .. .. .. .. .. .. ..@ nrows          : int 0
  .. .. .. .. .. .. ..@ elementType    : chr "ANY"
  .. .. .. .. .. .. ..@ elementMetadata: NULL
  .. .. .. .. .. .. ..@ metadata       : list()
  .. .. .. .. .. .. ..@ listData       :List of 2
  .. .. .. .. .. .. .. ..$ type       : chr(0) 
  .. .. .. .. .. .. .. ..$ description: chr(0) 
  .. .. .. .. ..@ metadata       : list()
  .. .. .. .. ..@ listData       : Named list()
  .. .. ..@ elementType    : chr "GRanges"
  .. .. ..@ metadata       : list()
  .. .. ..@ partitioning   :Formal class 'PartitioningByEnd' [package "IRanges"] with 5 slots
  .. .. .. .. ..@ end            : int [1:39376] 0 0 0 0 0 0 0 0 0 0 ...
  .. .. .. .. ..@ NAMES          : chr [1:39376] "100287102" "653635" "102466751" "107985730" ...
  .. .. .. .. ..@ elementType    : chr "ANY"
  .. .. .. .. ..@ elementMetadata: NULL
  .. .. .. .. ..@ metadata       : list()
  ..@ colData           :Formal class 'DFrame' [package "S4Vectors"] with 6 slots
  .. .. ..@ rownames       : chr [1:34] "GSM5073707" "GSM5073709" "GSM5073710" "GSM5073711" ...
  .. .. ..@ nrows          : int 34
  .. .. ..@ elementType    : chr "ANY"
  .. .. ..@ elementMetadata:Formal class 'DFrame' [package "S4Vectors"] with 6 slots
  .. .. .. .. ..@ rownames       : NULL
  .. .. .. .. ..@ nrows          : int 2
  .. .. .. .. ..@ elementType    : chr "ANY"
  .. .. .. .. ..@ elementMetadata: NULL
  .. .. .. .. ..@ metadata       : list()
  .. .. .. .. ..@ listData       :List of 2
  .. .. .. .. .. ..$ type       : chr [1:2] "input" "input"
  .. .. .. .. .. ..$ description: chr [1:2] "" ""
  .. .. ..@ metadata       : list()
  .. .. ..@ listData       :List of 2
  .. .. .. ..$ sample: chr [1:34] "GSM5073707" "GSM5073709" "GSM5073710" "GSM5073711" ...
  .. .. .. ..$ status: Factor w/ 2 levels "COVID-19 Negative",..: 1 1 1 1 2 2 2 2 2 2 ...
  ..@ assays            :Formal class 'SimpleAssays' [package "SummarizedExperiment"] with 1 slot
  .. .. ..@ data:Formal class 'SimpleList' [package "S4Vectors"] with 4 slots
  .. .. .. .. ..@ listData       :List of 1
  .. .. .. .. .. ..$ counts: int [1:39376, 1:34] 8 26 2 4 3 6 0 15 36 9 ...
  .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. .. .. .. ..$ : chr [1:39376] "100287102" "653635" "102466751" "107985730" ...
  .. .. .. .. .. .. .. ..$ : chr [1:34] "GSM5073707" "GSM5073709" "GSM5073710" "GSM5073711" ...
  .. .. .. .. ..@ elementType    : chr "ANY"
  .. .. .. .. ..@ elementMetadata: NULL
  .. .. .. .. ..@ metadata       : list()
  ..@ NAMES             : NULL
  ..@ elementMetadata   :Formal class 'DFrame' [package "S4Vectors"] with 6 slots
  .. .. ..@ rownames       : NULL
  .. .. ..@ nrows          : int 39376
  .. .. ..@ elementType    : chr "ANY"
  .. .. ..@ elementMetadata: NULL
  .. .. ..@ metadata       : list()
  .. .. ..@ listData       : Named list()
  ..@ metadata          :List of 1
  .. ..$ version:Classes 'package_version', 'numeric_version'  hidden list of 1
  .. .. ..$ : int [1:3] 1 44 0

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

Другой момент — переименование строк. Хотя в структуре объекта трудно понять наверняка, где расположены искомые строки, но нам это и не нужно. В общем случае достаточно использовать rownames(dds), чтобы их извлечь. Начальная матрица, у которой есть имена строк, никуда не делась (Картинка 2). В коде вы можете заметить, что мы не переименовываем строки, так как имена в таком численном формате entrez id, как на Картинке 2, нам еще пригодятся. Если же требуется лишь рассчитать дифференциальную экспрессию и построить визуализацию, то:

if (length(rownames(dds) == length(symbols))){
      rownames(dds) <- symbols
} else {
      print(»Длины векторов названий генов не совпадают»)
}

И смело используем код для визуализации после DESeq и results, например:

DESeq2::plotMA(dds, ylim = c(-5,5))

Каждая точка Картинке 4 — ген. Таких точек столько, сколько строк в матрице на Картинке 2. Положение точки по оси Y отражает, что происходит с геном в эксперименте по отношению к контролю - его экспрессия падает либо возрастает. Притом, здесь экспрессия это логарифмический, т.е. безразмерный показатель. Положение по оси X - это среднее количество транскриптов (экземпляров гена) в контроле. Таким способом мы можем быстро оценить, происходит ли в нашем эксперименте что-то, влияющее на экспрессию генов. Как видим, много синих точек - значимо дифференциально-экспрессирующихся генов.

Для более комплексных визуализаций — Heatmap:

# Извлечение нормализованных данных для значимых генов
significant_gene_names <- rownames(significant_genes)
normalized_counts <- counts(dds, normalized = TRUE)[significant_gene_names, ]

# Логарифмирование для уменьшения эффекта выбросов
log_normalized_counts <- log2(normalized_counts + 1)

# Создание аннотации образцов
sample_annotation <- HeatmapAnnotation(
  status = meta$status,
  col = list(status = c(«control» = «blue», «case» = «red»))
)

# Цветовая шкала для тепловой карты
col_fun <- colorRamp2(c(-2, 0, 2), c(«blue», «white», «red»))

# Построение тепловой карты
Heatmap(
  log_normalized_counts,
  name = «Expression»,
  top_annotation = sample_annotation,
  col = col_fun,
  show_row_names = TRUE,
  show_column_names = TRUE,
  row_names_gp = gpar(fontsize = 8),
  column_names_gp = gpar(fontsize = 8),
  cluster_rows = TRUE,
  cluster_columns = TRUE,
  row_title = «Significant Genes»,
  column_title = «Samples»,
  heatmap_legend_param = list(
    title = «Expression»,
    at = c(-2, 0, 2),
    labels = c(«Low», «Medium», «High»)
  )
)

ComplexHeatmap повзоляет строить визуализации дифференциальной экспрессии. При построении происходит и иерархическая кластеризация. Результат можно увидеть на Картинке 5. На ней гены (строки) сгруппированы по паттерну изменения экспрессии. Образцы (столбцы) тоже. Контроль — люди, не переболевшие коронавирусом - помечены над диаграммой сиреневыми прямоугольниками. Все остальные - переболевшие.

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

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

Другой способ визуализации - VolcanoPlot:

EnhancedVolcano(
  res,
  lab = rownames(res),                  # Метки генов
  x = 'log2FoldChange',                 # Ось X — log2FoldChange
  y = 'pvalue',                         # Ось Y — P-value (или можно использовать padj для скорректированных значений)
  pCutoff = 0.005,                       # Пороговое значение p-value
  FCcutoff = 1,                         # Порог log2FoldChange для отображения значительных генов
  title = 'Volcano plot of differentially expressed genes',
  xlim = c(-5, 5),                      # Ограничение для X (log2FoldChange)
  ylim = c(0, -log10(min(res$pvalue, na.rm = TRUE)))  # Ограничение для Y (p-value)
)

X — насколько значительно меняется экспрессия по модулю (сильно, слабо, уменьшается или увеличивается в эксперименте по сравнению с контролем), Y - статистическая значимость. Вертикальные отсечки - отсечки по модулю изменения экспрессии, который исследователь считает существенным, в данном случае 2 (шкала безразмерная). Горизонтальная отсечка - p-value = 0.005 по оси Y. Итого красные точки - статистически значимо и сильно изменившие свою экспрессию гены в о переболевших с симптомами нарушения обоняния по сравнению со здоровыми людьми, зеленым - сильно, но статистически незначимо, синим - статистически значимо, но очень слабо (статистически значимо изменений НЕТ), серым - ни то, ни другое

Volcano plot в данном случае более информативен, т.к. отражает как статистическую значимость, так и то, насколько изменилась экспрессия в ту или иную сторону, подсвечивает 30 наиболее значимо и значительно изменивших экспрессию генов.

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

Сигнальный путь — это каскад белков, который в наших клетках передает сигнал. Такой провод, где сигнал передается от белка к белку, благодаря чему могут возникать ответные реакции. Мы можем посмотреть, к каким сигнальным путям принадлежат найденные нами белки, и установить, где случается изменение. Для этого необходимы базы Reactome и KEGG. Именно поэтому выше было написано, что мы пока не будем менять наименования на общепринятые символы. База Reactome требует именно циферных обозначений генов, как на Картинке 2:

# Поскольку ваши гены уже в формате ENTREZ ID, используем rownames напрямую
entrez_ids <- rownames(significant_genes)
entrez_ids <- na.omit(entrez_ids)  # Убираем NA значения

# Выполняем обогащение по сигнальным и метаболическим путям с помощью KEGG
kegg_enrich <- enrichKEGG(gene = entrez_ids,
                          organism = 'hsa',   # hsa для человека
                          pvalueCutoff = 0.05)

# Анализ обогащения по путям Reactome
reactome_enrich <- enrichPathway(gene = entrez_ids,
                                 organism = «human»,
                                 pvalueCutoff = 0.05,
                                 readable = TRUE)  

# Визуализация результатов обогащения 

## Reactome pathways
dotplot(reactome_enrich, showCategory = 20) + ggtitle(«Reactome Pathway Enrichment»)

## KEGG pathways
dotplot(kegg_enrich, showCategory = 20) + ggtitle(«KEGG Pathway Enrichment»)

Картинка 7. Пример обогащения путей по Reactome.

Картинка 8. Пример обогащения путей по KEGG.

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

Далее, исходя из целей исследования, можно анализировать подобные диаграммы. Две вышеупомянутые картинки демонстрируют, что существенно изменяется экспрессия белков, связанных с передачей рецепторного сигнала. То есть происходит нарушение обоняния, так как анализировались данные обонятельного эпителия. База данных Reactome позволяет установить, какие конкретно сигнальные пути задействуются, либо быстро получить доступ к функциям конкретных белков, а также найти публикации по данной теме. Таким образом, обладая большим массивом данных секвенирования транскриптома, мы можем проанализировать, как клетки отвечают на внешние стимулы и посредством какие сигналов они это осуществляют. Исходя из дизайна эксперимента, можно проверять гипотезы по типу «Нарушение функции A связано с изменением экспрессии сигнальных путей процесса B внутри клетки, а не процесса C « и другие.

Заключение

Транскриптомный анализ bulk RNA-seq — относительно простая задача, однако требует базовых знаний в статистике, а также понимании синтаксиса R. Он прост еще в том, что, как правило, предобработанные матрицы экспрессии, которые можно непосредственно запустить в анализ, находятся в открытом доступе в Gene Expression Omnibus. В отличие от single-cell, порог входа в bulk RNA-seq гораздо ниже. На примере одного из датасетов, мы рассмотрели поэтапное проведение такого анализа. Конечно, можно выдвигать разные гипотезы для проверки, после чего код будет претерпевать изменения.

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

P.S. Статья написана по мотивам моего решения задачи Турнира Трех Наук. Отдельная благодарность Екатерине Олейник, члену оргкомитета Турнира, за активную помощь в работе над статьей и ее редактуре.

Транскриптомный анализ: как посчитать гены?