Язык R

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

Лекция 12

Анна Валяева

21 ноября 2025

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

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

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

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

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

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

Конец!