Язык R

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

Лекция 15

Анна Валяева

13 декабря 2024

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

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

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

  • Сборка транскриптома 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)

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

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

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

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

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

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

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

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

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

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

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

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")

DESeq

dds <- DESeq(dds)
  • считает нормировочные коэффициенты
  • для каждого гена строится 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)

Функциональная аннотация генов

Гадание по списку генов

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

 [1] "Pi15"          "Gm28438"       "Mettl21c"      "Mettl21e"     
 [5] "Col3a1"        "Zdbf2"         "2310015K22Rik" "Grem2"        
 [9] "Pfkfb3"        "Galnt5"        "A530058N18Rik" "Actc1"        
[13] "C130080G10Rik" "Sel1l2"        "Cdh4"          "B230312C02Rik"
[17] "Postn"         "Kcnab1"        "Mlf1"          "Col11a1"      
[21] "Abca4"         "Tmeff1"        "Mup10"         "Mup12"        
[25] "Efcab7"        "Pla2g5"        "Gm13091"       "Sema3a"       
[29] "Adamts3"       "Fras1"         "Mmd2"          "Col1a2"       
[33] "Pon1"          "Gm15737"       "Cacna2d4"      "2310014F06Rik"
[37] "Otoa"          "Adam12"        "Cyp2e1"        "Cilp2"        
[41] "Plet1os"       "AI118078"      "Tinag"         "Gm28979"      
[45] "Gm9795"        "Tet1"          "Aldh1l2"       "Kera"         
[49] "Wif1"          "Adcy1"         "2310065F04Rik" "Gm12300"      
[53] "Chad"          "Gm11752"       "Kcnf1"         "Rian"         
[57] "Irf4"          "Prss55"        "Pebp4"         "C1qtnf3"      
[61] "Gm17473"       "Sntb1"         "Wisp1"         "Rmi2"         
[65] "Tmem45a"       "Gm26885"       "H2-M9"         "Pla2g7"       
[69] "Nrep"          "Gm15337"       "Megf10"        "Aldh3b2"

Что делать дальше? Смотреть на каждый ген по отдельности?..

Список генов

Что про эти гены можем сказать?..

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

Хотим перенести анализ с уровня отдельных генов на уровень клеточных процессов.

Откуда брать аннотации для генов?

Можно создать самому:

TERM GENE
term1 gene1
term1 gene2
term2 gene2
term3 gene1
term3 gene3

GMT формат

covid19_gs <- clusterProfiler::read.gmt("COVID-19_gs.gmt")

https://software.broadinstitute.org/cancer/software/gsea/wiki/index.php/Data_formats

Много разных аннотаций доступны на сайте Enrichr.

Over-Representation Analysis

Анализ перепредставленности:

  • есть список генов
  • сопоставляем гены из списка и гены из категории
    • не оказалось ли в нашем списке таких генов больше, чем ожидалось бы по случайным причинам?
  • используем гипергеометрический тест (односторонний точный тест Фишера)


Список генов для ORA

Какой список подойдет?

  • дифф экспрессируемые гены
  • гены, значимо повысившие/понизившие экспрессию
  • любой другой список интересных/значимых генов
  • гены/белки/метаболиты, что угодно!

Гипергеометрический тест

\(p = 1 - \sum_{i=0}^{k-1}\frac{\binom{M}{i}\binom{N-M}{n-i}}{\binom{N}{n}}\)

  • N - всего генов (universe)
  • M - генов в категории
  • n - генов в списке
  • k - генов из категории в списке

Не забываем поправку на множественное тестирование!

Gene Set Enrichment Analysis

Анализ обогащения набора генов:

  • есть список всех [экспрессируемых] генов
  • список сортируем:
    • в начале - гены с повышенной экспрессией
    • в конце - гены с пониженной экспрессией
    • в середине - гены без изменения экспрессии
  • считаем скор - ES - обогащение генами из категории в начале или конце списка
  • рассчитываем p-значение с помощью пермутационного теста


Список генов для GSEA

Как ранжировать список?

  • по значениям log2 Fold Change
  • по комбинации p-значения (поправленного) и log2 Fold Change
    • например: \(| log_{10}(pval) | \space * \space sign(log_2(FoldChange))\)
  • …

Gene Set Variation Analysis

  • есть матрица экспрессии по всем образцам
  • внутри каждого образца ранжируем гены по уровню экспрессии
  • для каждого образца и набора генов рассчитываем скор - ES, сравнивая ранги генов, входящих в данный набор и вне него (K-S статистика)
  • получаем матрицу GSVA-скоров генных сигнатур (наборов генов)

Инструменты

Специальные пакеты

Онлайн-сервисы

clusterProfiler

clusterProfiler

BiocManager::install("clusterProfiler") 
library(clusterProfiler)

Книга-мануал по использованию пакета - Biomedical Knowledge Mining using GOSemSim and clusterProfiler

Пример данных

Результат сравнения экспрессии генов (анализ дифференциальной экспрессии генов) между молодыми и старыми мышами.

mouse
# A tibble: 15,615 × 3
   GeneName      padj log2FoldChange
   <chr>        <dbl>          <dbl>
 1 Rp1          0.957        0.0554 
 2 Sox17        0.994       -0.00283
 3 RP23-37D15.1 0.542        0.679  
 4 RP23-34E15.4 0.938       -0.0272 
 5 Mrpl15       0.288        0.319  
 6 RP23-34E15.5 0.782        0.256  
 7 Lypla1       0.943       -0.0372 
 8 RP24-426M1.3 0.385        0.192  
 9 Tcea1        0.705        0.180  
10 Gm6104       0.462        0.593  
# ℹ 15,605 more rows

Конвертация gene IDs

Разные пакеты и сервисы могут принимать только определенные идентификаторы генов. Не забывайте по возможное неоднозначное соответствие между разными идентификаторами.

В составе пакета clusterProfiler есть функция bitr():

# Аннотация для модельного организма - Mus musculus
BiocManager::install("org.Mm.eg.db") 

bitr(
  geneID = mouse$GeneName,           # вектор идентификаторов генов
  fromType = "SYMBOL",               # откуда
  toType = c("ENSEMBL", "ENTREZID"), # куда
  OrgDb = "org.Mm.eg.db")            # какой модельный организм
  SYMBOL            ENSEMBL ENTREZID
1    Rp1 ENSMUSG00000025900    19888
2  Sox17 ENSMUSG00000025902    20671
5 Mrpl15 ENSMUSG00000033845    27395

Какие варианты gene IDs доступны?

idType(OrgDb = "org.Mm.eg.db")
 [1] "ACCNUM"       "ALIAS"        "ENSEMBL"      "ENSEMBLPROT"  "ENSEMBLTRANS"
 [6] "ENTREZID"     "ENZYME"       "EVIDENCE"     "EVIDENCEALL"  "GENENAME"    
[11] "GENETYPE"     "GO"           "GOALL"        "IPI"          "MGI"         
[16] "ONTOLOGY"     "ONTOLOGYALL"  "PATH"         "PFAM"         "PMID"        
[21] "PROSITE"      "REFSEQ"       "SYMBOL"       "UNIPROT"

ORA с помощью clusterProfiler

Проведем анализ обогащения генов, которые значимо (поправленное p-значение ≤ 0.05) выше экспрессируются в молодых мышах по сравнению со старыми, используя аннотацию Gene Ontology, включающую все 3 ее части (BP + CC + MF).

mouse_up <- mouse %>% filter(padj <= 0.05, log2FoldChange > 0)

ora_up_go <- enrichGO(
  gene = mouse_up$GeneName,  # вектор идентификаторов генов
  OrgDb = "org.Mm.eg.db",    # модельный организм
  keyType = "SYMBOL",        # тип идентификатора генов
  ont = "ALL",               # аннотация GO BP+CC+MF
  universe = mouse$GeneName) # все детектированные гены

Результат ORA

Это объект S4 класса.

ora_up_go
#
# over-representation test
#
#...@organism    Mus musculus 
#...@ontology    GOALL 
#...@keytype     SYMBOL 
#...@gene    chr [1:241] "Adhfe1" "Pi15" "Gm28438" "Il1rl2" "Mettl21c" "Mettl21e" ...
#...pvalues adjusted by 'BH' with cutoff <0.05 
#...18 enriched terms found
'data.frame':   18 obs. of  10 variables:
 $ ONTOLOGY   : chr  "BP" "BP" "BP" "BP" ...
 $ ID         : chr  "GO:0030199" "GO:0030198" "GO:0043062" "GO:0045229" ...
 $ Description: chr  "collagen fibril organization" "extracellular matrix organization" "extracellular structure organization" "external encapsulating structure organization" ...
 $ GeneRatio  : chr  "8/189" "14/189" "14/189" "14/189" ...
 $ BgRatio    : chr  "51/12964" "243/12964" "243/12964" "243/12964" ...
 $ pvalue     : num  6.56e-07 1.29e-05 1.29e-05 1.29e-05 4.83e-09 ...
 $ p.adjust   : num  1.64e-03 8.08e-03 8.08e-03 8.08e-03 6.41e-07 ...
 $ qvalue     : num  1.63e-03 8.03e-03 8.03e-03 8.03e-03 6.14e-07 ...
 $ geneID     : chr  "Col3a1/Col5a2/Sfrp2/Col11a1/Col1a2/Adamts2/Loxl2/Adamts12" "Col3a1/Col5a2/Postn/Sfrp2/Col11a1/Col8a2/Adamts3/Col1a2/Antxr1/Adamts2/Loxl2/Adamts12/Sema5a/Kazald1" "Col3a1/Col5a2/Postn/Sfrp2/Col11a1/Col8a2/Adamts3/Col1a2/Antxr1/Adamts2/Loxl2/Adamts12/Sema5a/Kazald1" "Col3a1/Col5a2/Postn/Sfrp2/Col11a1/Col8a2/Adamts3/Col1a2/Antxr1/Adamts2/Loxl2/Adamts12/Sema5a/Kazald1" ...
 $ Count      : int  8 14 14 14 23 23 20 4 4 7 ...
#...Citation
 T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu.
 clusterProfiler 4.0: A universal enrichment tool for interpreting omics data.
 The Innovation. 2021, 2(3):100141

Результат ORA

Нам интереснее всего слот result.

# ora_up_go@result %>% head()
as.data.frame(ora_up_go) %>% head() 
           ONTOLOGY         ID                                   Description
GO:0030199       BP GO:0030199                  collagen fibril organization
GO:0030198       BP GO:0030198             extracellular matrix organization
GO:0043062       BP GO:0043062          extracellular structure organization
GO:0045229       BP GO:0045229 external encapsulating structure organization
GO:0031012       CC GO:0031012                          extracellular matrix
GO:0030312       CC GO:0030312              external encapsulating structure
           GeneRatio   BgRatio       pvalue     p.adjust       qvalue
GO:0030199     8/189  51/12964 6.555300e-07 1.644725e-03 1.634685e-03
GO:0030198    14/189 243/12964 1.287780e-05 8.077599e-03 8.028290e-03
GO:0043062    14/189 243/12964 1.287780e-05 8.077599e-03 8.028290e-03
GO:0045229    14/189 243/12964 1.287780e-05 8.077599e-03 8.028290e-03
GO:0031012    23/193 366/13035 4.832772e-09 6.413766e-07 6.135139e-07
GO:0030312    23/193 367/13035 5.090290e-09 6.413766e-07 6.135139e-07
                                                                                                                                                          geneID
GO:0030199                                                                                             Col3a1/Col5a2/Sfrp2/Col11a1/Col1a2/Adamts2/Loxl2/Adamts12
GO:0030198                                                  Col3a1/Col5a2/Postn/Sfrp2/Col11a1/Col8a2/Adamts3/Col1a2/Antxr1/Adamts2/Loxl2/Adamts12/Sema5a/Kazald1
GO:0043062                                                  Col3a1/Col5a2/Postn/Sfrp2/Col11a1/Col8a2/Adamts3/Col1a2/Antxr1/Adamts2/Loxl2/Adamts12/Sema5a/Kazald1
GO:0045229                                                  Col3a1/Col5a2/Postn/Sfrp2/Col11a1/Col8a2/Adamts3/Col1a2/Antxr1/Adamts2/Loxl2/Adamts12/Sema5a/Kazald1
GO:0031012 Col3a1/Col5a2/Col6a3/Angptl1/Itih2/Fbn1/Postn/Col11a1/Tnc/Col8a2/Adamts3/Fras1/Col1a2/Tinag/Igf1/Adamts2/Sparc/Chad/Thbs4/Loxl2/Gpc6/Adamts12/Kazald1
GO:0030312 Col3a1/Col5a2/Col6a3/Angptl1/Itih2/Fbn1/Postn/Col11a1/Tnc/Col8a2/Adamts3/Fras1/Col1a2/Tinag/Igf1/Adamts2/Sparc/Chad/Thbs4/Loxl2/Gpc6/Adamts12/Kazald1
           Count
GO:0030199     8
GO:0030198    14
GO:0043062    14
GO:0045229    14
GO:0031012    23
GO:0030312    23

Функции dplyr

Для работы с получаемыми в ходе анализа обогащений объектами в пакете clusterProfiler есть свои dplyr-функции: arrange(), filter(), mutate(), select(), slice(), summarise().

Они могут маскировать функции из пакета dplyr и выдавать ошибки при работе с датафреймами. Поэтому при вызове этих функций стоит указывать пакет: dplyr::filter().

ora_up_go %>% filter(ONTOLOGY == "BP") %>% head()
           ONTOLOGY         ID                                   Description
GO:0030199       BP GO:0030199                  collagen fibril organization
GO:0030198       BP GO:0030198             extracellular matrix organization
GO:0043062       BP GO:0043062          extracellular structure organization
GO:0045229       BP GO:0045229 external encapsulating structure organization
           GeneRatio   BgRatio      pvalue    p.adjust      qvalue
GO:0030199     8/189  51/12964 6.55530e-07 0.001644725 0.001634685
GO:0030198    14/189 243/12964 1.28778e-05 0.008077599 0.008028290
GO:0043062    14/189 243/12964 1.28778e-05 0.008077599 0.008028290
GO:0045229    14/189 243/12964 1.28778e-05 0.008077599 0.008028290
                                                                                                         geneID
GO:0030199                                            Col3a1/Col5a2/Sfrp2/Col11a1/Col1a2/Adamts2/Loxl2/Adamts12
GO:0030198 Col3a1/Col5a2/Postn/Sfrp2/Col11a1/Col8a2/Adamts3/Col1a2/Antxr1/Adamts2/Loxl2/Adamts12/Sema5a/Kazald1
GO:0043062 Col3a1/Col5a2/Postn/Sfrp2/Col11a1/Col8a2/Adamts3/Col1a2/Antxr1/Adamts2/Loxl2/Adamts12/Sema5a/Kazald1
GO:0045229 Col3a1/Col5a2/Postn/Sfrp2/Col11a1/Col8a2/Adamts3/Col1a2/Antxr1/Adamts2/Loxl2/Adamts12/Sema5a/Kazald1
           Count
GO:0030199     8
GO:0030198    14
GO:0043062    14
GO:0045229    14

Функции для ORA

  • enrichGO() - анализ обогащения Gene Ontology категорий
  • enrichKEGG() - анализ KEGG путей
  • enrichWP() - анализ WikiPathways путей
  • ReactomePA::enrichPathway() - анализ Reactome путей
  • enrichDO() - анализ Disease Ontology, Network of Cancer Gene и Disease Gene Network
  • enrichMeSH() - анализ MeSH категорий
  • enricher() - анализ любой кастомной аннотации

GSEA с помощью clusterProfiler

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

Для работы с KEGG нам нужны идентификаторы генов ENTREZ.

mouse_ids <- bitr(
  geneID = mouse$GeneName, 
  fromType = "SYMBOL",  # из символьных обозначений
  toType = "ENTREZID",  # получаем Entrez IDs
  OrgDb = "org.Mm.eg.db") %>% 
  # Сразу объединяем с таблицей с log2 Fold change и p-значениями
  inner_join(mouse, by = c("SYMBOL" = "GeneName"))

GSEA с помощью clusterProfiler

# Создаем отсортированный список генов
gene_list <- mouse_ids$log2FoldChange
names(gene_list) <- mouse_ids$ENTREZID
gene_list <- sort(gene_list, decreasing = TRUE)

gsea_kegg <- gseKEGG(
  geneList = gene_list, # список генов
  organism = "mmu",     # модельный организм
  pvalueCutoff = 0.1)   # фильтр по p-значению

Результат GSEA

gsea_kegg
#
# Gene Set Enrichment Analysis
#
#...@organism    mmu 
#...@setType     KEGG 
#...@geneList    Named num [1:14337] 4.01 3.84 3.31 2.95 2.76 ...
 - attr(*, "names")= chr [1:14337] "75104" "676527" "100303644" "100039054" ...
#...nPerm    
#...pvalues adjusted by 'BH' with cutoff <0.1 
#...6 enriched terms found
'data.frame':   6 obs. of  11 variables:
 $ ID             : chr  "mmu05169" "mmu05150" "mmu04060" "mmu04061" ...
 $ Description    : chr  "Epstein-Barr virus infection - Mus musculus (house mouse)" "Staphylococcus aureus infection - Mus musculus (house mouse)" "Cytokine-cytokine receptor interaction - Mus musculus (house mouse)" "Viral protein interaction with cytokine and cytokine receptor - Mus musculus (house mouse)" ...
 $ setSize        : int  171 38 144 47 71 71
 $ enrichmentScore: num  -0.448 -0.627 -0.452 -0.606 0.435 ...
 $ NES            : num  -1.82 -1.99 -1.8 -2.01 1.75 ...
 $ pvalue         : num  6.31e-05 3.83e-04 3.61e-04 9.80e-04 1.29e-03 ...
 $ p.adjust       : num  0.0209 0.0423 0.0423 0.0811 0.0851 ...
 $ qvalue         : num  0.0187 0.0378 0.0378 0.0724 0.0761 ...
 $ rank           : num  3201 1038 2566 2737 1733 ...
 $ leading_edge   : chr  "tags=38%, list=22%, signal=30%" "tags=32%, list=7%, signal=29%" "tags=31%, list=18%, signal=26%" "tags=49%, list=19%, signal=40%" ...
 $ core_enrichment: chr  "57296/12367/14991/14964/29857/15024/17246/18709/19184/19094/12567/12018/18707/22030/21354/22031/15015/22059/239"| __truncated__ "12260/12262/14961/16408/12268/14960/12259/14129/20345/15894/14969/16668" "230979/12161/12774/22163/21937/16194/12977/13051/85030/100504346/242316/20301/12156/18829/12772/16154/14560/148"| __truncated__ "12978/230979/12774/21937/16194/12977/13051/100504346/20301/18829/12772/16154/14825/76527/17000/57266/56066/5598"| __truncated__ ...
#...Citation
 T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu.
 clusterProfiler 4.0: A universal enrichment tool for interpreting omics data.
 The Innovation. 2021, 2(3):100141

setReadable()

Конвертация Entrez IDs в символьные обозначения в столбце core_enrichment:

gsea_kegg <- setReadable(gsea_kegg, OrgDb = "org.Mm.eg.db", keyType = "ENTREZID")

gsea_kegg %>% select(ID, core_enrichment) %>% head()
               ID
mmu05169 mmu05169
mmu05150 mmu05150
mmu04060 mmu04060
mmu04061 mmu04061
mmu04974 mmu04974
mmu04512 mmu04512
                                                                                                                                                                                                                                                                                                                                                                                         core_enrichment
mmu05169 Psmd8/Casp3/H2-M3/H2-D1/Mapk12/H2-T10/Mdm2/Pik3r2/Psmc5/Mapk11/Cdk4/Bak1/Pik3cd/Traf2/Tap1/Traf3/H2-Q4/Trp53/Oas1b/Irf3/Psmc4/H2-DMa/Ccne1/Akt2/H2-Q6/Oas2/Irf7/Tyk2/H2-DMb2/Jun/Rela/Tradd/Psmc3/Jak3/Psmd4/Psmd3/Adrm1/H2-T22/Ncor2/Tab1/Nfkbib/H2-Ab1/Itgal/Isg15/H2-Q7/Bax/H2-Aa/Oas1a/Ccnd1/Akt1/Bid/Relb/Rb1/Icam1/H2-Eb1/Tlr2/Map3k14/Nfkbie/Cdkn1a/E2f2/Tnf/Nfkb2/Gadd45a/H2-T25/Cxcl10
mmu05150                                                                                                                                                                                                                                                                                                                           C1qb/C1qc/H2-Ab1/Itgal/C4b/H2-Aa/C1qa/Fcgr1/Selplg/Icam1/H2-Eb1/Krt18
mmu04060                                                                                                          Tnfrsf14/Bmp6/Ccr5/Tnfrsf4/Tnfrsf1a/Il6ra/Csf1/Cx3cr1/Tnfrsf25/Gm13304/Gdf6/Ccl27a/Bmp2/Ccl21a/Ccr2/Il10ra/Gdf10/Cxcl1/Edar/Il34/Thpo/Il3ra/Tgfb3/Relt/Ltbr/Tnfrsf13b/Cxcl14/Gdf9/Inhbb/Csf3r/Cxcl11/Amhr2/Cxcr6/Ngfr/Cntfr/Il17re/Cxcl13/Ccl8/Tnf/Il20rb/Cx3cl1/Lep/Eda2r/Cxcl10/Gdf5
mmu04061                                                                                                                                                                                                                                             Csf1r/Tnfrsf14/Ccr5/Tnfrsf1a/Il6ra/Csf1/Cx3cr1/Gm13304/Ccl27a/Ccl21a/Ccr2/Il10ra/Cxcl1/Il34/Ltbr/Cxcl14/Cxcl11/Cxcl13/Ccl8/Tnf/Il20rb/Cx3cl1/Cxcl10
mmu04974                                                                                                                                                                                                          Col11a1/Col1a2/Col3a1/Slc36a2/Ace2/Col12a1/Col1a1/Col6a5/Col9a1/Col8a2/Col5a2/Col26a1/Kcne3/Col5a1/Col6a3/Col6a6/Slc36a4/Col6a1/Xpnpep2/Col4a3/Col4a4/Col5a3/Dpp4/Kcnj13/Atp1a4/Col8a1
mmu04512                                                                                                                                                                                                                                                                                      Col1a2/Fras1/Chad/Frem1/Col1a1/Col6a5/Col9a1/Tnc/Npnt/Thbs4/Lama3/Itga6/Col6a3/Col6a6/Col6a1/Col4a3/Col4a4

Функции для GSEA

  • gseGO() - анализ обогащения Gene Ontology категорий
  • gseKEGG() - анализ KEGG путей
  • gseWP() - анализ WikiPathways путей
  • ReactomePA::gsePathway() - анализ Reactome путей
  • gseDO(), gseNCG(), gseDGN() - анализ Disease Ontology, Network of Cancer Gene и Disease Gene Network
  • gseMeSH() - анализ MeSH категорий
  • GSEA() - анализ любой кастомной аннотации

Визуализация обогащений

Столбчатая диаграмма

barplot(ora_up_go) +
  labs(title = "Activated \n GO categories")

Дотплот

dotplot(ora_down_go) +
  labs(title = "Suppressed \n GO categories")

Граф генов в обогащенных категориях

cnetplot(ora_up_go)

Граф генов в обогащенных категориях

cnetplot(ora_up_go, foldChange = gene_list)

Граф обогащенных категорий

ora_down_go %>% filter(ONTOLOGY == "BP") %>% 
  enrichplot::pairwise_termsim() %>%  # сходство между категориями
  emapplot()

Дерево обогащенных GO категорий

ora_down_go %>% filter(ONTOLOGY == "MF") %>% 
  enrichplot::pairwise_termsim() %>%  # сходство между категориями
  enrichplot::treeplot(nCluster= 3)

Ridgeline plot

ridgeplot(gsea_kegg)

GSEA running score

gseaplot(gsea_kegg, geneSetID = 1, title = gsea_kegg$Description[1])

GSEA running score

enrichplot::gseaplot2(gsea_kegg, geneSetID = 2, title = gsea_kegg$Description[2])

pathview

BiocManager::install("pathview")
library("pathview")

hsa04110 <- pathview(
  # Вектор log2 fold change с именами генов 
  gene.data = gene_list, 
  pathway.id = "mmu05150",
  species = "mmu",
  limit = list(gene = c(-1,1), cpd = 1),
  low = list(gene = "blue", cpd = "green"))

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

Про анализ обогащений

Конец!