Enrichment analysis

Material

Exercises

tum_vs_norm <- readRDS("day3/tum_vs_norm_day3-2.rds")
limma_de <- readRDS("day3/limma_de_day3-2.rds")
proB <- readRDS("day3/proB_day3-2.rds")

Load the following packages:

If the FindMarkers or FindAllMarkers functions were used, we obtained a table listing only the significant genes, but we don’t have any information of fold change for the non-significant genes. Therefore, we can use the over-representation analysis which is a threshold-based method. Using our list of significant genes, we can test if any gene set is over-represented among significant genes or not using a test similar to a Fisher test to compare differences in proportions.

The clusterProfiler package provides functions for over-representation analysis of Gene Ontology gene sets (among other functions, including functions for actual GSEA) or KEGG gene sets.

Genes can be labeled using different types of labels, eg symbol, Ensembl ID, Entrez ID. To list the allowed label types use:

# BiocManager::install("org.Hs.eg.db", update = FALSE)
library(org.Hs.eg.db)
AnnotationDbi::keytypes(org.Hs.eg.db)

About OrgDb

For other organisms, you can find available OrgDbs at bioconductor

Let’s select a set of genes that are downregulated in the tumor cells compared to normal:

tum_down  <- subset(limma_de,
                    limma_de$logFC < -1 
                      & limma_de$adj.P.Val <  0.05)
tum_down_genes <- rownames(tum_down)

We can do a Gene Ontology term over-representation analysis based on this set of genes. Make sure you check out the help of this function to understand its arguments:

?enrichGO
tum_vs_norm_go <- clusterProfiler::enrichGO(tum_down_genes,
                                            "org.Hs.eg.db",
                                            keyType = "SYMBOL",
                                            ont = "BP",
                                            minGSSize = 50)

The results are stored in the @result slot:

View(tum_vs_norm_go@result)
ID Description GeneRatio BgRatio RichFactor FoldEnrichment oddsRatio
GO:0007059 GO:0007059 chromosome segregation 96/805 420/18842 0.2285714 5.349991 7.402393
GO:0000070 GO:0000070 mitotic sister chromatid segregation 64/805 205/18842 0.3121951 7.307305 10.962223
GO:0000819 GO:0000819 sister chromatid segregation 69/805 246/18842 0.2804878 6.565157 9.459746
GO:0098813 GO:0098813 nuclear chromosome segregation 76/805 315/18842 0.2412698 5.647213 7.763532
GO:0140014 GO:0140014 mitotic nuclear division 72/805 281/18842 0.2562278 5.997321 8.378859
GO:0000280 GO:0000280 nuclear division 82/805 422/18842 0.1943128 4.548126 5.903320
NoteThe columns GeneRatio and BgRatio

The columns GeneRatio and BgRatio that are in the enrichResult object represent the numbers that are used as input for the Fisher’s exact test.

The two numbers (M/N) in the GeneRatio column are:

  • M: Number of genes of interest (in our case tum_down_genes) that are in the GO set
  • N: Number of genes of interest with any GO annoation.

The two numbers (k/n) in the BgRatio column are:

  • k: Number of genes in the universe that are in the GO set
  • n: Number of genes in the universe with any GO annoation

A low p-value resulting from the Fisher’s exact means that M/N is signficantly greater than k/n. 

Some GO terms seem redundant because they contain many of the same genes, which is a characteristic of Gene Ontology gene sets. We can simplify this list by removing redundant gene sets:

enr_go <- clusterProfiler::simplify(tum_vs_norm_go)
View(enr_go@result)
ID Description GeneRatio BgRatio RichFactor FoldEnrichment oddsRatio
GO:0007059 GO:0007059 chromosome segregation 96/805 420/18842 0.2285714 5.349991 7.402393
GO:0000070 GO:0000070 mitotic sister chromatid segregation 64/805 205/18842 0.3121951 7.307305 10.962223
GO:0098813 GO:0098813 nuclear chromosome segregation 76/805 315/18842 0.2412698 5.647213 7.763532
GO:0000280 GO:0000280 nuclear division 82/805 422/18842 0.1943128 4.548126 5.903320
GO:0010948 GO:0010948 negative regulation of cell cycle process 61/805 309/18842 0.1974110 4.620644 5.881075
GO:0044772 GO:0044772 mitotic cell cycle phase transition 71/805 434/18842 0.1635945 3.829127 4.709670

We can quite easily generate a plot called an enrichment map with the enrichplot package:

enrichplot::emapplot(enrichplot::pairwise_termsim(enr_go),
                     showCategory = 30)

Instead of testing for Gene Ontology terms, we can also test for other gene set collections. For example the Hallmark collection from MSigDB:

library(msigdbr)
# library(msigdbdf)
gmt <- msigdbr::msigdbr(species = "human", category = "H")
Warning: The `category` argument of `msigdbr()` is deprecated as of msigdbr 10.0.0.
ℹ Please use the `collection` argument instead.

We can use the function enricher to test for over-representation of any set of genes of the Hallmark collection. We have to include the “universe”, i.e. the full list of background, non significant genes, against which to test for differences in proportions:

tum_vs_norm_enrich <- clusterProfiler::enricher(gene = tum_down_genes,
                                                universe = rownames(proB),
                                                pAdjustMethod = "BH",
                                                pvalueCutoff  = 0.05,
                                                qvalueCutoff  = 0.05,
                                                TERM2GENE = gmt[,c("gs_name", "gene_symbol")])

When using the genes down-regulated in tumor, among the over-represented Hallmark gene sets, we have HALLMARK_G2M_CHECKPOINT, which includes genes involved in the G2/M checkpoint in the progression through the cell division cycle.

View(tum_vs_norm_enrich@result[tum_vs_norm_enrich@result$p.adjust < 0.05,])
ID Description GeneRatio BgRatio RichFactor FoldEnrichment oddsRatio
HALLMARK_E2F_TARGETS HALLMARK_E2F_TARGETS HALLMARK_E2F_TARGETS 80/346 195/3855 0.4102564 4.570920 8.876103
HALLMARK_G2M_CHECKPOINT HALLMARK_G2M_CHECKPOINT HALLMARK_G2M_CHECKPOINT 67/346 187/3855 0.3582888 3.991917 6.782049
HALLMARK_MITOTIC_SPINDLE HALLMARK_MITOTIC_SPINDLE HALLMARK_MITOTIC_SPINDLE 47/346 196/3855 0.2397959 2.671715 3.544702
HALLMARK_MYC_TARGETS_V1 HALLMARK_MYC_TARGETS_V1 HALLMARK_MYC_TARGETS_V1 34/346 193/3855 0.1761658 1.962772 2.296001
HALLMARK_ESTROGEN_RESPONSE_LATE HALLMARK_ESTROGEN_RESPONSE_LATE HALLMARK_ESTROGEN_RESPONSE_LATE 28/346 164/3855 0.1707317 1.902227 2.183777

Clear your environment:

gc()