Here's one 'on the house', in part copied from 2 previous answers that I have given on this topic:

rankedGenes <- c(1:20)
names(rankedGenes) <- c("BRCA1","BRCA2","TP53", "ATM", "ERCC1",
  "BRCC3", "TP63", "ERCC2", "ERCC3", "ERCC4", "TP63",
  "S100A2", "CXCL2", "CXCL3", "MS4A1", "MS4A2", "VCAM1",
  "VSX1", "CD8", "S100A2")

rankedGenes
 BRCA1  BRCA2   TP53    ATM  ERCC1  BRCC3   TP63  ERCC2  ERCC3  ERCC4   TP63 
     1      2      3      4      5      6      7      8      9     10     11 
S100A2  CXCL2  CXCL3  MS4A1  MS4A2  VCAM1   VSX1    CD8 S100A2 
    12     13     14     15     16     17     18     19     20

So, in my genes, BRCA1 is ranked highest... S100A2 is lowest.

require(topGO)
require(org.Hs.eg.db)

selection <- function(x) TRUE
allGO2genes <- annFUN.org(
  whichOnto = "BP",
  feasibleGenes = NULL,
  mapping = "org.Hs.eg.db",
  ID = "symbol")
GOdata <- new(
  "topGOdata",
  ontology = "BP", # also available are 'CC' and 'MF'
  allGenes = rankedGenes,
  annot = annFUN.GO2genes,
  GO2genes = allGO2genes,
  geneSel = selection,
  nodeSize = 5)

Here, we make use of the rank information by using the Kolmogorov-Smirnov (K-S) test on the enriched terms

results.ks <- runTest(GOdata, algorithm = "classic", statistic = "ks")
goEnrichment <- GenTable(
  GOdata,
  KS = results.ks,
  orderBy = "KS",
  topNodes = 100,
  numChar = 99)
goEnrichment$KS <- as.numeric(goEnrichment$KS)
goEnrichment <- goEnrichment[goEnrichment$KS < 0.05,] # filter terms for KS p<0.05
goEnrichment <- goEnrichment[,c("GO.ID","Term","KS")]
goEnrichment
require(ggplot2)
library(scales)

ntop <- 30
ggdata <- goEnrichment[1:ntop,]
ggdata$Term <- factor(ggdata$Term, levels = rev(ggdata$Term)) # fixes order
gg1 <- ggplot(ggdata,
  aes(x = Term, y = -log10(KS), size = -log10(KS), fill = -log10(KS))) +

  expand_limits(y = 1) +
  geom_point(shape = 21) +
  scale_size(range = c(2.5,12.5)) +
  scale_fill_continuous(low = 'royalblue', high="red4") +

  xlab('') + ylab('Enrichment score') +
  labs(
    title="GO Biological processes",
    subtitle="Top 30 terms ordered by Kolmogorov-Smirnov p-value",
    caption = 'Cut-off lines drawn at equivalents of p=0.05, p=0.01, p=0.001') +

  geom_hline(yintercept = c(-log10(0.05), -log10(0.01), -log10(0.001)),
    linetype = c("dotted", "longdash", "solid"),
    colour = c("black", "black", "black"),
    size = c(0.5, 1.5, 3)) +

  theme_bw(base_size = 24) +
  theme(
    legend.position = 'right',
    legend.background = element_rect(),
    plot.title = element_text(angle = 0, size = 16, face="bold", vjust = 1),
    plot.subtitle = element_text(angle = 0, size = 14, face="bold", vjust = 1),
    plot.caption = element_text(angle = 0, size = 12, face="bold", vjust = 1),

    axis.text.x = element_text(angle = 0, size = 12, face="bold", hjust = 1.10),
    axis.text.y = element_text(angle = 0, size = 12, face="bold", vjust = 0.5),
    axis.title = element_text(size = 12, face="bold"),
    axis.title.x = element_text(size = 12, face="bold"),
    axis.title.y = element_text(size = 12, face="bold"),
    axis.line = element_line(colour="black"),

    #Legend
    legend.key = element_blank(), # removes the border
    legend.key.size = unit(1, "cm"), # Sets overall area/size of the legend
    legend.text = element_text(size = 14, face = "bold"), # Text size
    title = element_text(size = 14, face = "bold")) +

  coord_flip()

hhh



Source link