Язык R

и его применение в биоинформатике

Лекция 12

Анна Валяева

24 ноября 2025

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

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

Пример задачи: как изучаемое вещество X влияет на профиль экспрессии генов в клеточной линии HeLa?

На какие вопросы поможет ответить анализ дифф экспрессии генов?

  • химическое или физическое воздействие vs контроль
  • больные vs здоровые
  • разные органы, ткани, клеточные типы
  • этапы развития
  • деплеция гена vs контроль
  • оверэкспрессия гена vs контроль
  • …

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

Kukurba and Montgomery, 2015

О чем задуматься перед РНК-секвенированием?

  • poly(A)-селекция или rRNA деплеция
  • одно- или парно-концевые риды [single-end / paired-end]
  • цепь-специфичная библиотека или нет
  • количество биологических реплик
  • глубина секвенирования (30 млн чтений по рекомендациям ENCODE)
  • минимизация батч эффекта

Experimental planning considerations

Batch effect

Hicks et al., 2015 (bioRxiv)

Как из прочтений получить экспрессию?

Как из прочтений получить экспрессию?

Для работы “с нуля” необходимы:

  • прочтения (fastq)
  • геном или транскриптом
  • разметка генов (gtf/gff)
  • информация о дизайне эксперимента (что есть что)

Gene expression matrix


Дизайн эксперимента

Задача: узнать, экспрессия каких генов изменилась в ответ на воздействие (treatment) по сравнению с контролем (control).

При этом в каждой группе как минимум 2 биологических повторности!

Как проверить, что ген X в среднем имеет экспрессию в одной группе образцов выше или ниже, чем в другой?

Сравнивать между образцами можно только нормированную экспрессию.

Нормализация

  • размер библиотеки (глубина секвенирования)

Нормализация

  • размер библиотеки (глубина секвенирования)
  • состав транскриптома (ДЭ гены, количество экспрессируемых генов, контаминации)

Нормализация

  • размер библиотеки (глубина секвенирования)
  • состав транскриптома
  • длина гена

Варианты нормализации

  • размер библиотеки (глубина секвенирования)

CPM (Counts Per Million)

\(CPM = \frac{Number \space of \space reads \space mapped \space to \space gene \space × \space 10^6}{Total \space number \space of \space mapped \space reads \space}\)

  • размер библиотеки (глубина секвенирования) + длина гена

RPKM (Reads Per Kilobase Million), FPKM (Fragments Per Kilobase Million)

\(RPKM = \frac{Number \space of \space reads \space mapped \space to \space gene \space × 10^3 × \space 10^6}{Total \space number \space of \space mapped \space reads \space × Gene \space length \space in \space bp}\)

TPM (Transcripts Per Kilobase Million)

\(TPM = A \space × \space \frac{1}{\sum(A)} \space × \space 10^6\), где \(A = \frac{Total \space reads \space mapped \space to \space gene \space × \space 10^3}{Gene \space length \space in \space bp}\)

Как проверить, что ген X в среднем имеет экспрессию в одной группе образцов выше или ниже, чем в другой?

  • Средний уровень экспрессии в группе CTRL = 5

  • Средний уровень экспрессии в группе TRT = 10

  • \(Fold \space Change = \frac{mean \space expession_{treatment}} {mean \space expession_{control}} = \frac{10} {5} = 2\)

  • \(log_2FoldChange = log_2(\frac{mean \space expession_{treatment}} {mean \space expession_{control}}) = log_2(\frac{10} {5}) = 1\)


Экспрессия повышена в 2 раза: \(log_2FoldChange = 1\)

Экспрессия не изменилась: \(log_2FoldChange = 0\)

Экспрессия понижена в 2 раза: \(log_2FoldChange = -1\)

Как проверить, что ген X в среднем имеет экспрессию в одной группе образцов выше или ниже, чем в другой?

\(log_2FoldChange = 1\)


Но насколько статистически значим этот результат?

Не могли ли мы его наблюдать по случайным причинам, например, из-за технического или биологического шума?

Статистическая значимость

Проведем статистический тест Манна-Уитни. В результате теста получили p-значение 0.056.

Что это значит?


P-значение - вероятность получить результат такой же, как видим сейчас, или еще экстремальней, при условии что нулевая гипотеза верна.

Гипотезы

Нулевая гипотеза

\(H_0\) - принимаемое по умолчанию скептическое предположение о том, что эффекта нет.

В контексте анализа дифф экспрессии: \(H_0\) - настоящее значение \(log_2FoldChange\) равно 0.

Альтернативная гипотеза

\(H_1\), или \(H_A\) - альтернативное предположение о том, что эффект есть.

О каком результате идет речь?

  • о выборке
  • о статистике, посчитанной по выборке:
    • среднее
    • медиана
    • \(log_2FoldChange\)
    • t-статистика
    • …

Более экстремальный результат - более маловероятный.

Уровень значимости

Если p-значение меньше заданного уровня значимости α (обычно 0.05), то нулевая гипотеза отвергается в пользу альтернативной - эффект есть.

Если p-значение больше α, то нет оснований отвергнуть нулевую гипотезу.


Уровень значимости выбирается заранее. Он контролирует вероятность ошибки I рода - вероятность ошибочно отвергнуть истинную нулевую гипотезу.

Проблема множественного тестирования

Наша задача - найти дифф экспрессируемые гены среди всех ненулевых генов.

Всего ненулевых генов ~10,000. К каждому гену мы применяем статистический тест. Даже если никаких настоящих изменений в экспрессии всех генов нет, просто из-за случайных причин примерно 5% тестов дадут p < 0.05.

Проблема множественного тестирования: чем больше тестов - тем больше шанс случайно найти “значимый” результат.

Нужно применять поправки на множественное тестирование!

Специальные подходы для анализа RNA-seq

Для анализа дифф экспрессии по данным РНК секвенирования разработаны специальные алгоритмы, учитывающие:

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

Если выборки большие, например, большие когорты больных и здоровых людей или десятки тысяч клеток в эксперименте по секвенированию единичных клеток, можно использовать универсальные тесты, например, критерий Манна-Уитни (Wilcoxon rank-sum test).

R пакеты для анализа дифф экспрессии

  • DESeq2
  • edgeR
  • limma + voom
  • sleuth
  • ballgown
  • NOISeq
  • DEXSeq
  • …

DESeq2

Bioconductor - DESeq2

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("DESeq2")

library(DESeq2)

DESeq2 workflow

dds <- DESeqDataSetFromMatrix(
  countData = counts_mtx,
  colData = coldata,
  design= ~ batch + condition)

dds <- DESeq(dds) # вся магия

res <- results(dds, name = "condition_trt_vs_untrt")

Матрица каунтов

counts <- read_tsv("../data/mouse_genes/GSE184348_counts.tsv")
str(counts, give.attr = FALSE)
spc_tbl_ [43,346 × 10] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
 $ GeneId  : chr [1:43346] "ENSMUSG00000102693" "ENSMUSG00000064842" "ENSMUSG00000051951" "ENSMUSG00000102851" ...
 $ GeneName: chr [1:43346] "RP23-271O17.1" "Gm26206" "Xkr4" "RP23-317L18.1" ...
 $ Young_1 : num [1:43346] 0 0 2 0 0 0 0 4 0 0 ...
 $ Young_2 : num [1:43346] 0 0 3 0 0 0 0 1 0 1 ...
 $ Young_3 : num [1:43346] 0 0 1 0 0 0 0 3 0 0 ...
 $ Young_4 : num [1:43346] 0 0 2 0 0 0 0 0 0 0 ...
 $ Old_1   : num [1:43346] 0 0 2 0 0 0 0 1 0 1 ...
 $ Old_2   : num [1:43346] 0 0 0 0 0 0 0 0 0 0 ...
 $ Old_3   : num [1:43346] 0 0 2 0 0 0 0 0 0 0 ...
 $ Old_4   : num [1:43346] 0 0 1 0 0 0 0 0 0 0 ...

Преобразовать в матрицу

counts_mtx <- counts %>% select(where(is.numeric)) %>% as.matrix()
rownames(counts_mtx) <- counts$GeneName

Метаданные

Создать датафрейм вручную или прочитать из файла

coldata <- data.frame(
  Mouse = str_remove(colnames(counts_mtx), "_\\d"),
  row.names = colnames(counts_mtx))
coldata
        Mouse
Young_1 Young
Young_2 Young
Young_3 Young
Young_4 Young
Old_1     Old
Old_2     Old
Old_3     Old
Old_4     Old

DESeqDataSet

  • матрица каунтов
  • датафрейм с метаданными (названия строк = названия образцов в матрице каунтов)
  • формула (design) для модели, включающая переменные из метаданных: ~ some_variable + batch + variable_of_interest
dds <- DESeqDataSetFromMatrix(
  countData = counts_mtx,
  colData = coldata,
  design = ~ Mouse)

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

tximport

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

# Создаем таблицу соотвествия транскрипт-ген
library(TxDb.Mmusculus.UCSC.mm10.knownGene)
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
k <- keys(txdb, keytype = "TXNAME")
tx2gene <- select(txdb, k, "GENEID", "TXNAME")

# Получаем список файлов kallisto
k_files <- list.files(
  file.path(data_dir, "kallisto"), 
  pattern = "abundance.tsv", 
  recursive = TRUE, full.names = TRUE)
names(k_files) <- paste0("sample", 1:6)

# Агрегируем экспрессию транскриптов
BiocManager::install("tximport")
library(tximport)
txi <- tximport(k_files, type = "kallisto", tx2gene = tx2gene)

# Создаем DESeq-объект
dds <- DESeqDataSetFromTximport(
  txi, 
  colData = coldata,
  design = ~ Mouse)

Нормализация DESeq2

\[s_{j} = \underset{i}{median}{k_{ij} \over (\prod_{v=1}^m k_{iv})^{1/m}}\]

  1. Среднее геометрическое для каждого гена (устойчиво к выбросам) - псевдо-референс
  2. Отношение экспрессии к псевдо-референсу
  3. Медиана отношений по каждому образцу

Основано на допущении: лишь относительно небольшое число генов действительно дифференциально экспрессируются.

Нормировочные коэффициенты

\(Normalized\ counts_{j} = {Raw\ counts_{j} \over s_{j}}\)

dds <- estimateSizeFactors(dds)
sizeFactors(dds) %>% as.numeric()
[1] 1.0488580 1.1423301 0.6885799 0.8511961 1.3068077 1.1323248 0.9218329
[8] 1.0793817
counts(dds, normalized = TRUE)[1:10, 1:2] 
               Young_1   Young_2
RP23-271O17.1 0.000000 0.0000000
Gm26206       0.000000 0.0000000
Xkr4          1.906836 2.6262111
RP23-317L18.1 0.000000 0.0000000
RP23-317L18.4 0.000000 0.0000000
RP23-317L18.3 0.000000 0.0000000
RP23-115I1.6  0.000000 0.0000000
Gm1992        3.813671 0.8754037
RP23-115I1.5  0.000000 0.0000000
RP23-115I1.2  0.000000 0.8754037

Успешная нормализация


Нет ли среди образцов выбросов?

  • корреляция
  • кластеризация
  • PCA

Трансформация каунтов

Для визуализации каунты трансформируют, чтобы сделать дисперсию “примерно одинаковой” для всех генов, независимо от их уровня экспрессии:

  • rlog() - regularized log transformation
  • vst() - variance stabilizing transformation

PCA в DESeq2

transformed_dds <- rlog(dds)
# или 
transformed_dds <- vst(dds)

plotPCA(transformed_dds, intgroup=c("condition"))

PCA

Сильно отличающиеся образцы - выбросы

Поэтому лучше секвенировать больше повторностей на всякий случай.

Batch effect in PCA

Переменную batch можно добавить в формулу для модели DESeq2.

Kucukural et al., 2019

PCA

plotPCA(vst(dds), intgroup = "Mouse")

Оценка дисперсии

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

dds <- estimateDispersions(dds)
plotDispEsts(dds)

DESeq

dds <- DESeq(dds)
  • считает нормировочные коэффициенты (size factors)
  • для каждого гена строится GLM
  • каунты моделируются с помощью отрицательного биномиального распределения: \(K_{ij} \sim NB(mean = \mu_{ij}, dispersion = \alpha_i)\)
  • дисперсия оценивается сначала для каждого гена, затем корректируется с помощью эмпирического байесовского подхода
  • подсчитываются log2 fold change значения
  • значимость изменений экспрессии оценивается с помощью теста Вальда (Wald) и проводится поправка на множественное тестирование

DESeq результаты

resultsNames(dds)
[1] "Intercept"          "Mouse_Young_vs_Old"
res <- results(dds, name = "Mouse_Young_vs_Old")

DESeq результаты

str(res)
Formal class 'DESeqResults' [package "DESeq2"] with 7 slots
  ..@ priorInfo      : list()
  ..@ rownames       : chr [1:43346] "RP23-271O17.1" "Gm26206" "Xkr4" "RP23-317L18.1" ...
  ..@ nrows          : int 43346
  ..@ elementType    : chr "ANY"
  ..@ elementMetadata:Formal class 'DFrame' [package "S4Vectors"] with 6 slots
  .. .. ..@ rownames       : NULL
  .. .. ..@ nrows          : int 6
  .. .. ..@ elementType    : chr "ANY"
  .. .. ..@ elementMetadata: NULL
  .. .. ..@ metadata       : list()
  .. .. ..@ listData       :List of 2
  .. .. .. ..$ type       : chr [1:6] "intermediate" "results" "results" "results" ...
  .. .. .. ..$ description: chr [1:6] "mean of normalized counts for all samples" "log2 fold change (MLE): Mouse Young vs Old" "standard error: Mouse Young vs Old" "Wald statistic: Mouse Young vs Old" ...
  ..@ metadata       :List of 6
  .. ..$ filterThreshold: Named num 4.45
  .. .. ..- attr(*, "names")= chr "60.27064%"
  .. ..$ filterTheta    : num 0.603
  .. ..$ filterNumRej   :'data.frame':  50 obs. of  2 variables:
  .. .. ..$ theta : num [1:50] 0.401 0.412 0.423 0.435 0.446 ...
  .. .. ..$ numRej: num [1:50] 733 744 750 754 755 768 785 795 798 807 ...
  .. ..$ lo.fit         :List of 2
  .. .. ..$ x: num [1:50] 0.401 0.412 0.423 0.435 0.446 ...
  .. .. ..$ y: num [1:50] 734 741 748 756 764 ...
  .. ..$ alpha          : num 0.1
  .. ..$ lfcThreshold   : num 0
  ..@ listData       :List of 6
  .. ..$ baseMean      : num [1:43346] 0 0 1.62 0 0 ...
  .. ..$ log2FoldChange: num [1:43346] NA NA 0.895 NA NA ...
  .. ..$ lfcSE         : num [1:43346] NA NA 1.15 NA NA ...
  .. ..$ stat          : num [1:43346] NA NA 0.78 NA NA ...
  .. ..$ pvalue        : num [1:43346] NA NA 0.435 NA NA ...
  .. ..$ padj          : num [1:43346] NA NA NA NA NA NA NA NA NA NA ...

DESeq результаты

res %>% as.data.frame() %>% 
  mutate(gene = rownames(.)) %>% 
  select(gene, log2FoldChange, padj) %>% 
  drop_na(padj) %>% 
  filter(padj < 0.05) %>% 
  arrange(desc(abs(log2FoldChange))) %>% head()
           gene log2FoldChange         padj
Krt8       Krt8      -7.854023 1.173674e-06
Krt18     Krt18      -7.791874 6.866451e-23
Astl       Astl      -5.908037 3.180374e-04
Fgf21     Fgf21      -5.892606 1.383096e-05
Prrt3     Prrt3      -5.576513 1.103255e-02
Gm15935 Gm15935      -5.530609 2.381310e-05

MA-plot

Volcano plot

EnhancedVolcano: publication-ready volcano plots with enhanced colouring and labeling

Что почитать?

Про DESeq2

Про анализ дифф. экспрессии

Данные: GSE184348 (Zhang et al., 2022)

Конец!