Disease Gene Set Enrichment Analysis
2017-10-30
Contents
1 GSEA algorithm
A common approach in analyzing gene expression profiles was identifying differential expressed genes that are deemed interesting. The enrichment analysis we demonstrated in Disease enrichment analysis vignette were based on these differential expressed genes. This approach will find genes where the difference is large, but it will not detect a situation where the difference is small, but evidenced in coordinated way in a set of related genes. Gene Set Enrichment Analysis (GSEA)1 directly addresses this limitation. All genes can be used in GSEA; GSEA aggregates the per gene statistics across genes within a gene set, therefore making it possible to detect situations where all genes in a predefined set change in a small but coordinated way. Since it is likely that many relevant phenotypic differences are manifested by small but consistent changes in a set of genes.
Genes are ranked based on their phenotypes. Given a priori defined set of gens S (e.g., genes shareing the same DO category), the goal of GSEA is to determine whether the members of S are randomly distributed throughout the ranked gene list (L) or primarily found at the top or bottom.
There are three key elements of the GSEA method:
- Calculation of an Enrichment Score.
- The enrichment score (ES) represent the degree to which a set S is over-represented at the top or bottom of the ranked list L. The score is calculated by walking down the list L, increasing a running-sum statistic when we encounter a gene in S and decreasing when it is not. The magnitude of the increment depends on the gene statistics (e.g., correlation of the gene with phenotype). The ES is the maximum deviation from zero encountered in the random walk; it corresponds to a weighted Kolmogorov-Smirnov-like statistic1.
- Esimation of Significance Level of ES.
- The p-value of the ES is calculated using permutation test. Specifically, we permute the gene labels of the gene list L and recompute the ES of the gene set for the permutated data, which generate a null distribution for the ES. The p-value of the observed ES is then calculated relative to this null distribution.
- Adjustment for Multiple Hypothesis Testing.
- When the entire gene sets were evaluated, DOSE adjust the estimated significance level to account for multiple hypothesis testing and also q-values were calculated for FDR control.
We implemented GSEA algorithm proposed by Subramanian1. Alexey Sergushichev implemented an algorithm for fast GSEA analysis in the fgsea2 package.
In DOSE3, user can use GSEA algorithm implemented in DOSE or fgsea by specifying the parameter by="DOSE" or by="fgsea". By default, DOSE use fgsea since it is much more fast.
1.1 Leading edge analysis and core enriched genes
Leading edge analysis reports Tags to indicate the percentage of genes contributing to the enrichment score, List to indicate where in the list the enrichment score is attained and Signal for enrichment signal strength.
It would also be very interesting to get the core enriched genes that contribute to the enrichment.
DOSE supports leading edge analysis and report core enriched genes in GSEA analysis.
1.2 gseDO fuction
In the following example, in order to speedup the compilation of this document, only gene sets with size above 120 were tested and only 100 permutations were performed.
library(DOSE)
data(geneList)
y <- gseDO(geneList,
nPerm = 100,
minGSSize = 120,
pvalueCutoff = 0.2,
pAdjustMethod = "BH",
verbose = FALSE)
head(y, 3)## ID Description setSize enrichmentScore
## DOID:0060084 DOID:0060084 cell type benign neoplasm 439 -0.2837045
## DOID:5614 DOID:5614 eye disease 450 -0.3125247
## DOID:5679 DOID:5679 retinal disease 299 -0.3676313
## NES pvalue p.adjust qvalues rank
## DOID:0060084 -1.254662 0.01250000 0.06922399 0.02923977 1890
## DOID:5614 -1.380387 0.01315789 0.06922399 0.02923977 1768
## DOID:5679 -1.591370 0.01315789 0.06922399 0.02923977 1768
## leading_edge
## DOID:0060084 tags=23%, list=15%, signal=20%
## DOID:5614 tags=22%, list=14%, signal=19%
## DOID:5679 tags=24%, list=14%, signal=21%
## core_enrichment
## DOID:0060084 37/268/2735/6752/3082/5914/2878/7508/100507436/3791/1301/1543/1027/1028/1958/7173/7201/8743/7450/596/947/2034/11197/4314/324/3964/196/3675/595/6843/573/2246/3912/1902/7507/2308/8829/6817/2934/1277/3953/55384/8773/8991/2057/4311/2247/6414/4588/5243/5468/358/1012/6469/9457/1733/81029/3952/126/79068/7058/4488/1009/4313/56034/3625/2944/6925/2099/3480/6387/5159/57496/1462/1289/2690/6571/4487/1191/54361/10752/5744/4922/2487/247/7043/1811/2952/367/3572/1287/4582/7031/3479/6424/4629/652/10647/5241/4969
## DOID:5614 3082/5914/2878/4153/3791/23247/1543/80184/6750/1958/2098/7450/596/9187/2034/482/948/1490/1280/3931/5737/4314/4881/2261/3426/187/629/6403/7042/6785/7507/2934/5176/4060/1277/7078/5950/2057/727/10516/4311/2247/1295/358/10203/2192/582/10218/57125/3485/585/1675/6310/2202/4313/2944/4254/3075/1501/2099/3480/4653/6387/3305/1471/857/4016/1909/4053/6678/1296/7033/4915/55812/1191/5654/10631/2152/2697/7043/2952/6935/2200/3572/7177/7031/3479/2006/10451/9370/771/3117/125/652/4693/5346/1524
## DOID:5679 2878/3791/23247/80184/6750/7450/596/9187/2034/482/948/1490/1280/5737/4314/4881/3426/187/629/6403/6785/2934/5176/7078/5950/727/10516/4311/2247/1295/358/10203/582/10218/57125/585/1675/6310/2202/4313/2944/4254/3075/2099/3480/4653/6387/1471/857/4016/1909/4053/6678/1296/4915/55812/1191/5654/10631/2697/2952/6935/2200/3479/2006/10451/9370/771/652/4693/5346/15241.3 gseNCG fuction
ncg <- gseNCG(geneList,
nPerm = 100,
minGSSize = 120,
pvalueCutoff = 0.2,
pAdjustMethod = "BH",
verbose = FALSE)
ncg <- setReadable(ncg, 'org.Hs.eg.db')
head(ncg, 3)## ID Description setSize enrichmentScore NES pvalue
## breast breast breast 133 -0.4869070 -1.873013 0.01388889
## lung lung lung 173 -0.3880662 -1.532406 0.01408451
## lymphoma lymphoma lymphoma 188 0.2999589 1.240661 0.03571429
## p.adjust qvalues rank leading_edge
## breast 0.04225352 0.02223870 2930 tags=33%, list=23%, signal=26%
## lung 0.04225352 0.02223870 2775 tags=31%, list=22%, signal=25%
## lymphoma 0.07142857 0.03759398 2087 tags=21%, list=17%, signal=18%
## core_enrichment
## breast KMT2A/ERBB3/SETD2/ARID1A/GPS2/NCOR1/RB1/MAP2K4/NF1/TP53/PIK3R1/STK11/CDKN1B/PTGFR/APC/CCND1/TRAF5/MAP3K1/ESR1/TBX3/FOXA1/GATA3
## lung SETD2/ATXN3L/LRP1B/BRD3/ARID1A/INHBA/RB1/ADCY1/LYRM9/NF1/CTNNB1/TP53/SATB2/STK11/CTIF/CTNNA3/KDR/COL11A1/FLT3/APC/ADGRL3/FGFR3/NCAM2/DIP2C/APLNR/SLIT2/EPHA3/RUNX1T1/ZMYND10/ZFHX4/GLI3/TNN/PLSCR4/DACH1/ERBB4
## lymphoma DUSP2/EZH2/PRDM1/MYC/ZWILCH/IKZF3/PLCG2/IDH2/HIST1H1C/MAGEC3/CD79B/ETV6/HIST1H1E/HIST1H1B/IRF8/CD28/SLC29A2/DUSP9/TNFAIP3/DNMT3A/SYK/TNF/BCR/HIST1H1D/DSC3/UBE2A/PABPC11.4 gseDGN fuction
dgn <- gseDGN(geneList,
nPerm = 100,
minGSSize = 120,
pvalueCutoff = 0.2,
pAdjustMethod = "BH",
verbose = FALSE)
dgn <- setReadable(dgn, 'org.Hs.eg.db')
head(dgn, 3)## ID Description setSize
## umls:C0029456 umls:C0029456 Osteoporosis 375
## umls:C0338656 umls:C0338656 Impaired cognition 342
## umls:C1272641 umls:C1272641 Systemic arterial pressure 318
## enrichmentScore NES pvalue p.adjust qvalues
## umls:C0029456 -0.3439046 -1.517108 0.01388889 0.09310987 0.05423242
## umls:C0338656 -0.3266625 -1.450775 0.01408451 0.09310987 0.05423242
## umls:C1272641 -0.3277208 -1.444783 0.01408451 0.09310987 0.05423242
## rank leading_edge
## umls:C0029456 1766 tags=23%, list=14%, signal=20%
## umls:C0338656 1997 tags=23%, list=16%, signal=20%
## umls:C1272641 1758 tags=23%, list=14%, signal=20%
## core_enrichment
## umls:C0029456 HGF/PTH1R/CYP1A1/JAG1/ROR2/FLT3/CUL9/EEF1A2/THSD4/BCL2/ITGAV/WIF1/GREM2/COL15A1/HPGDS/VGLL3/SLIT3/NRIP1/TMEM135/MGP/PLCL1/OSBPL1A/PIBF1/SELP/SPRY1/MMP13/ID4/SPP2/COL1A2/AOX1/ARHGEF3/GSN/TSC22D3/ATP1B1/NR5A2/ANKH/COL1A1/LEPR/THSD7A/GC/FGF2/PPARG/NOX4/ZNF266/GHRH/BHLHE40/SLC19A2/THBD/FLNB/KL/LEP/HSD17B4/CTSK/FTO/MMP2/ESR1/IGF1R/PTN/IRAK3/HSPA1L/CST3/GHR/SPARC/KDM4B/LRP1/INPP4B/BMPR1B/PTHLH/DPT/FRZB/GSTT1/AR/TNFRSF11B/IRS1/WLS/GSTM3/TGFBR3/TPH1/IGF1/SFRP4/CORIN/BMP4/CHAD/FOXA1/PGR
## umls:C0338656 NR3C1/CAPN3/SLC2A10/CREBBP/ZNF224/ITM2B/ELK3/CLN5/GAD1/BACE1/HGF/SERPINA3/MBL2/SST/EGR1/INSR/UTRN/ARL4D/PVALB/EEF1A2/DYM/CD36/RAB40AL/RBMS3/TREM2/PER3/OXTR/TSC1/CDR1/IGFALS/TPPP/SELP/NGF/BCHE/KCNS3/APBB2/TRPM4/RUNX1T1/MME/ABCB1/PPARG/MVP/NME8/SPG11/LPL/SLC26A4/FHL5/KL/LEP/FTO/NAIP/SORL1/ESR1/ABCC8/CST3/LAMA2/HHAT/LRP1/CLU/ALB/SPON1/NTS/HTRA1/GSTT1/GRIA2/MAGI2/IRS1/TAT/COL4A5/AASS/IGF1/ITPR1/BMP4/LRP2/MAPT/ERBB4/GRP
## umls:C1272641 HGF/SERPINA3/SCNN1A/HECTD4/MBL2/CYP1A1/ENOX1/STK39/INSR/IGFBP7/TENM4/ITGAV/CRBN/FHIT/PRKG1/CD36/CTGF/GULP1/SLC24A3/NPY2R/MMP3/AHR/LONRF1/FBXL7/MEF2C/RPS6KA2/ID4/CDADC1/TGFBR2/VTCN1/CYP7A1/CITED2/SCNN1G/FAM155A/NEFL/JAM2/RGS17/PPARG/CDH13/SYNE1/GRK4/LPL/TRPS1/SEMA5A/KL/UCP1/HSD17B4/KCNMA1/PRLR/CDH11/SORL1/SGCD/NOV/ESR1/IGF1R/ZNF652/ACBD4/FGF18/FRY/ZNF385D/AGTR1/TNFRSF11B/IRS1/DCN/CARTPT/ADRA2A/TBX3/ADIPOQ/CORIN/ADH1B/NOVA1/LRP22 Visualization
2.1 cnetplot
cnetplot(ncg, categorySize="pvalue", foldChange=geneList)
2.2 enrichMap
enrichMap(y, n=20)
2.3 gseaplot
gseaplot(y, geneSetID = y$ID[1], title=y$Description[1])
References
1. Subramanian, A. et al. Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proceedings of the National Academy of Sciences of the United States of America 102, 15545–15550 (2005).
2. S., A. An algorithm for fast preranked gene set enrichment analysis using cumulative statistic calculation. biorxiv doi:10.1101/060012
3. Yu, G., Wang, L.-G., Yan, G.-R. & He, Q.-Y. DOSE: An r/bioconductor package for disease ontology semantic and enrichment analysis. Bioinformatics 31, 608–609 (2015).