Язык R

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

Лекция 10

Анна Валяева

7 ноября 2025

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

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

Зачем нужен РНК-сек?

  • Сборка транскриптома de novo
  • сравнение уровней экспрессии
    • больные 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 биологических повторности!

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

  • Тест на то, что разница в средней экспресии гена между группами равна 0

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

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

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)

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

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

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

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

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

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

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

\[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

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


Другие варианты нормализации

  • глубина секвенирования
  • длина гена

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}\)

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

  • корреляция
  • кластеризация
  • 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)

Конец!