32 Functional Enrichment (GSEA)
Gene set enrichment analysis (fgsea) tests whether a ranked list of genes
is non-randomly distributed within curated gene sets, using a running-sum
statistic rather than a hard significance cutoff — genes just below a
p-value threshold still contribute to the result. msigdbr provides the
MSigDB gene set collections in a ready-to-use table; the C7 (immunologic
signatures) collection is used here as it is the most directly relevant to
a PASC vs. Healthy immune comparison.
32.2 Ranked gene list
fgsea expects a named numeric vector (gene symbol -> rank statistic),
built with tibble::deframe() from a two-column data frame. The original
example used a small placeholder list of random ranks; the ranked list here
instead comes from a real result already computed in the book — the
Wilcoxon per-cell-type regulon test’s avg_log2FC, restricted to the RNA
marker genes underlying cluster_0’s regulons (Regulon Differential
Expression per Cell Type chapter also computes de_result for the same
cluster, if that chapter’s object is still in memory; it is recomputed here
so this chapter can run independently).
library(presto)
DefaultAssay(combined) <- "RNA"
seurat_object_celltype <- combined[, combined$cell_labels == "cluster_0"]
Idents(seurat_object_celltype) <- seurat_object_celltype$individual_condition
de_result_rna <- FindMarkers(seurat_object_celltype,
ident.1 = "Healthy",
ident.2 = "PASC",
logfc.threshold = 0,
min.pct = 0)
ranks <- de_result_rna %>%
rownames_to_column("genes") %>%
dplyr::select(genes, ranks = avg_log2FC) %>%
deframe()32.3 Run fgsea against MSigDB C7
scoreType = "std" is used rather than the original "pos", since
avg_log2FC ranks are signed (up in Healthy vs. up in PASC) — "pos"
assumes a non-negative statistic such as a rank of p-values.
fgseaRes <- fgsea(fgsea_sets, stats = ranks, scoreType = "std")
fgseaRes %>%
arrange(padj) %>%
head(20)
#> pathway
#> <char>
#> 1: NAKAYA_PBMC_FLUMIST_AGE_18_50YO_3DY_DN
#> 2: GSE9988_ANTI_TREM1_VS_LOW_LPS_MONOCYTE_UP
#> 3: GSE9988_ANTI_TREM1_VS_LPS_MONOCYTE_UP
#> 4: GSE14769_UNSTIM_VS_60MIN_LPS_BMDM_DN
#> 5: GSE9988_ANTI_TREM1_VS_CTRL_TREATED_MONOCYTES_UP
#> 6: GSE9988_LOW_LPS_VS_ANTI_TREM1_AND_LPS_MONOCYTE_DN
#> 7: GSE9988_LPS_VS_LPS_AND_ANTI_TREM1_MONOCYTE_DN
#> 8: HARALAMBIEVA_PBMC_FLUARIX_AGE_50_74YO_CORR_WITH_28D_MEM_B_CELL_RESPONSE_AT_28DY_POSITIVE
#> 9: NAKAYA_PBMC_FLUMIST_AGE_18_50YO_7DY_DN
#> 10: GSE7460_CTRL_VS_TGFB_TREATED_ACT_FOXP3_MUT_TCONV_DN
#> 11: GSE46606_UNSTIM_VS_CD40L_IL2_IL5_DAY1_STIMULATED_BCELL_DN
#> 12: GSE7460_FOXP3_MUT_VS_WT_ACT_TCONV_UP
#> 13: ANDERSON_BLOOD_CN54GP140_ADJUVANTED_WITH_GLA_AF_AGE_18_45YO_1DY_DN
#> 14: ANDERSON_BLOOD_CN54GP140_ADJUVANTED_WITH_GLA_AF_AGE_18_45YO_1DY_UP
#> 15: ANDERSON_BLOOD_CN54GP140_ADJUVANTED_WITH_GLA_AF_AGE_18_45YO_3DY_DN
#> 16: ANDERSON_BLOOD_CN54GP140_ADJUVANTED_WITH_GLA_AF_AGE_18_45YO_3DY_UP
#> 17: ANDERSON_BLOOD_CN54GP140_ADJUVANTED_WITH_GLA_AF_AGE_18_45YO_6HR_DN
#> 18: ANDERSON_BLOOD_CN54GP140_ADJUVANTED_WITH_GLA_AF_AGE_18_45YO_6HR_UP
#> 19: ANDERSON_BLOOD_CN54GP140_ADJUVANTED_WITH_GLA_AF_AGE_18_45YO_7DY_DN
#> 20: ANDERSON_BLOOD_CN54GP140_ADJUVANTED_WITH_GLA_AF_AGE_18_45YO_7DY_UP
#> pathway
#> <char>
#> pval padj log2err ES NES
#> <num> <num> <num> <num> <num>
#> 1: 2.566624e-05 0.1339264 0.57561026 -0.4925682 -1.3797669
#> 2: 5.222759e-05 0.1362618 0.55733224 -0.6206277 -1.5851544
#> 3: 9.821683e-05 0.1708318 0.53843410 -0.6038104 -1.5422187
#> 4: 8.266619e-04 0.5073423 0.47727082 -0.5646105 -1.4417019
#> 5: 7.478516e-04 0.5073423 0.47727082 -0.5675901 -1.4497069
#> 6: 8.631768e-04 0.5073423 0.47727082 -0.5722452 -1.4578650
#> 7: 4.821455e-04 0.5073423 0.49849311 -0.5895436 -1.4991387
#> 8: 8.750634e-04 0.5073423 0.47727082 -0.4279628 -1.2334445
#> 9: 5.834471e-04 0.5073423 0.47727082 -0.4484921 -1.2720696
#> 10: 1.260182e-03 0.6575630 0.45505987 -0.5522460 -1.4083958
#> 11: 1.698088e-03 0.8055113 0.45505987 -0.5617289 -1.4346754
#> 12: 2.064472e-03 0.8977011 0.43170770 -0.5628867 -1.4363680
#> 13: 3.041894e-01 1.0000000 0.10208011 -0.4596798 -1.0678446
#> 14: 6.129032e-01 1.0000000 0.07235709 0.3984168 0.9223240
#> 15: 9.046653e-01 1.0000000 0.05090829 -0.3712462 -0.7463068
#> 16: 9.127789e-01 1.0000000 0.05049830 -0.3671477 -0.7380677
#> 17: 2.313663e-02 1.0000000 0.35248786 -0.6152808 -1.3724335
#> 18: 8.683694e-01 1.0000000 0.05132233 0.3809787 0.7853204
#> 19: 5.776398e-01 1.0000000 0.07343814 0.7000534 0.9929392
#> 20: 4.347826e-01 1.0000000 0.08603472 0.4828065 1.0081286
#> pval padj log2err ES NES
#> <num> <num> <num> <num> <num>
#> size
#> <int>
#> 1: 663
#> 2: 192
#> 3: 193
#> 4: 199
#> 5: 193
#> 6: 191
#> 7: 190
#> 8: 1202
#> 9: 825
#> 10: 195
#> 11: 194
#> 12: 198
#> 13: 83
#> 14: 76
#> 15: 33
#> 16: 33
#> 17: 61
#> 18: 33
#> 19: 7
#> 20: 36
#> size
#> <int>
#> leadingEdge
#> <list>
#> 1: FAM210B,PLK3,QSOX1,LMNA,VEGFA,DDIT4,...[319]
#> 2: TBC1D2,ZNF668,STX3,TRIB1,SLC37A2,PPARG,...[71]
#> 3: TBC1D2,MAPKAPK3,STX3,TRIB1,SLC37A2,PPARG,...[85]
#> 4: HILPDA,PHLDA1,OSM,ZC3H12C,ATF3,F3,...[105]
#> 5: TBC1D2,STX3,SRXN1,TRIB1,DDIT4,PHLDA1,...[66]
#> 6: TBC1D2,STX3,TRIB1,AGAP3,PPARG,KLHL21,...[68]
#> 7: TBC1D2,STX3,TRIB1,AGAP3,PPARG,KLHL21,...[66]
#> 8: C15orf39,TBC1D2,SLC8A1,ST6GALNAC2,TTYH3,CXCL16,...[457]
#> 9: SOCS3,TRIB1,PHLDA1,RFLNB,IFITM1,KLHL21,...[385]
#> 10: DAPK2,TTC39B,NRARP,GGN,GPR68,CTSW,...[64]
#> 11: ZMYND19,MCM10,ETV5,SLFN11,HILPDA,PAQR7,...[70]
#> 12: TOX2,CYP11A1,ASB2,XKR8,STK16,UCP2,...[78]
#> 13: DPEP2,ARID3A,ARAP1,STX2,ARF3,DVL3,...[43]
#> 14: NRSN2,STAR,SYP,CRB3,FSTL4,SIGLEC6,...[20]
#> 15: KCTD2,DPEP2,ACTG1,DNAJB12,RAB7A,ATP6V1A,...[9]
#> 16: BMP8B,GCK,EML2
#> 17: DPF1,DDIT4,CMTM3,SLC39A13,SGCA,SLC35G6,...[21]
#> 18: SLC4A4,AP1M2,PTGR2,RASSF6,RAPSN,UFSP2,...[9]
#> 19: UQCR10,NUDT16L1,PRKCZ
#> 20: STAR,SLC26A4-AS1,SIGLEC6,RASSF6,ZNF568,TMEM156,...[9]
#> leadingEdge
#> <list>32.4 Run fgsea against a custom gene set list
fgsea also accepts any user-defined named list of gene sets in the same
format as fgsea_sets above. The regulon-to-target-gene table
(TF_regulation.txt, Regulon Activity Scoring chapter) is used here as a
custom set: each regulon’s target genes form one gene set, letting the same
enrichment test ask whether a regulon’s targets are shifted as a group,
independent of the MSigDB collections.
TF_regulation <- read.table("data/TF_regulation.txt", header = TRUE)
regulon_sets <- split(TF_regulation$item, TF_regulation$name)
fgseaRes_regulons <- fgsea(regulon_sets, stats = ranks, scoreType = "std")
fgseaRes_regulons %>%
arrange(padj) %>%
head(20)
#> pathway pval padj log2err ES
#> <char> <num> <num> <num> <num>
#> 1: ARID3A(+) 0.22265625 0.9300691 0.1275053 -0.8477407
#> 2: BCL11A(-) 0.15517241 0.9300691 0.1541910 -0.8974939
#> 3: CCDC88A(+) 0.04755548 0.9300691 0.3217759 -0.6885660
#> 4: CREB3L2(+) 0.24110672 0.9300691 0.1226792 -0.7094758
#> 5: EBF1(-) 0.15134100 0.9300691 0.1563124 -0.8153245
#> 6: EGR1(+) 0.24952015 0.9300691 0.1182875 -0.7471150
#> 7: ELF1(+) 0.18181818 0.9300691 0.1437590 -0.7416574
#> 8: ETS1(+) 0.19354839 0.9300691 0.1357409 -0.9006584
#> 9: ETS1(-) 0.12701613 0.9300691 0.1766943 0.9313188
#> 10: ETV5(+) 0.02126481 0.9300691 0.3524879 0.9825996
#> 11: FOSL1(+) 0.12723658 0.9300691 0.1752040 -0.7678558
#> 12: FOSL2(+) 0.12079208 0.9300691 0.1797823 -0.6224690
#> 13: HDX(+) 0.08893281 0.9300691 0.2114002 -0.9526597
#> 14: HIC1(+) 0.25616698 0.9300691 0.1157344 -0.8784475
#> 15: HIVEP1(+) 0.27075099 0.9300691 0.1147507 -0.8889742
#> 16: HLX(+) 0.26538462 0.9300691 0.1142665 -0.8210443
#> 17: IRF7(+) 0.25494071 0.9300691 0.1188150 -0.8943984
#> 18: IRF9(+) 0.27470356 0.9300691 0.1137873 -0.8870977
#> 19: JDP2(+) 0.16666667 0.9300691 0.1492075 -0.7786425
#> 20: JUNB(+) 0.18379447 0.9300691 0.1429011 -0.9121544
#> pathway pval padj log2err ES
#> <char> <num> <num> <num> <num>
#> NES size
#> <num> <int>
#> 1: -1.209724 7
#> 2: -1.234467 5
#> 3: -1.408207 37
#> 4: -1.189298 15
#> 5: -1.290430 11
#> 6: -1.203686 12
#> 7: -1.243244 15
#> 8: -1.208857 4
#> 9: 1.267068 3
#> 10: 1.336836 3
#> 11: -1.300812 16
#> 12: -1.258002 35
#> 13: -1.269846 3
#> 14: -1.179046 4
#> 15: -1.184957 3
#> 16: -1.206542 8
#> 17: -1.192187 3
#> 18: -1.182456 3
#> 19: -1.275870 13
#> 20: -1.215855 3
#> NES size
#> <num> <int>
#> leadingEdge
#> <list>
#> 1: SLC8A1,SOX5,BCL11B,PIK3R5,CD86
#> 2: PPARG,PID1,ARHGAP26
#> 3: SLC8A1,IGF2BP2,PPARG,RFX2,ANKDD1A,ANXA1,...[26]
#> 4: SLC8A1,PPARG,PID1,DMXL2,EHD4,MAML3,...[9]
#> 5: WDFY3,LCP2,ITGA5,PIK3R5,LRP1
#> 6: SLC8A1,DOCK4,DMXL2,ITGA5,PLCB1,ZEB2
#> 7: PPARG,PID1,DOCK4,UBASH3B,ANTXR2,GNAQ,...[8]
#> 8: SLCO3A1
#> 9: RAPH1,CTBP2
#> 10: SVIL
#> 11: PPARG,PID1,LCP2,ITGA5,EHD4,PLCB1,...[10]
#> 12: PPARG,RFX2,LCP2,ITGA5,ANKDD1A,ANXA1,...[13]
#> 13: SLC8A1
#> 14: ITGA5,EHD4,TIMP2
#> 15: C1orf21
#> 16: PPARG,GAS7,CD86,SRGN,GNAQ,AHR
#> 17: FGD4,CD86
#> 18: PID1,LRP1,SRGN
#> 19: PPARG,LCP2,MAML3,GAS7,ZMIZ1,LRP1,...[7]
#> 20: BCL11B
#> leadingEdge
#> <list>