-
Notifications
You must be signed in to change notification settings - Fork 0
/
GS_Pathway_analysis.Rmd
59 lines (41 loc) · 1.34 KB
/
GS_Pathway_analysis.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
---
title: "GS + Pathway enrichment"
author: "AEKS"
date: "2023-11-08"
output: html_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
```{r packages}
library(DESeq2)
library(tidyverse)
library(fgsea)
library(qusage)
library(enrichR)
```
```{r GSEA}
#rank gene list
ranked_genes <- -log10(DESeq2_Output$FDR)*DESeq2_Output$logFC
names(ranked_genes) <- DESeq2_Output$gene_name
#import chosen gene set
pathway_gene_sets <- read.gmt("h.all.v7.5.1.symbols.gmt")
fgsea_results <- fgseaMultilevel(pathways = pathway_gene_sets,
stats = ranked_genes)
#plot gene sets
ggplot(fgsea_results, aes(reorder(pathway, NES), NES)) +
geom_col(aes(fill = padj < 0.05)) +
coord_flip() +
labs(x = "Pathway", y = "Normalized Enrichment Score (NES)")
```
```{r pathway enrichment}
#check database list
dbsList <- listEnrichrDbs()
#choose pathway databases
dbs <- c("Reactome_2022", "KEGG_2021_Human", "BioPlanet_2019", "GO_Biological_Process_2023", "GO_Cellular_Component_2023", "GO_Molecular_Function_2023", "MSigDB_Hallmark_2020", "GWAS_Catalog_2023")
#Filter FDR
DESeq2_Output_filtered <- DESeq2_Output[DESeq2_Output$FDR < 0.05, ]
enriched_pathways <- enrichr(DESeq2_Output_filtered$gene_name, dbs)
#Extract gene names
enriched_pathways_gene_names <- enriched_pathways[["Pathway name"]][c(), c(), drop = FALSE]
```