rGREAT.md

Pathway analysis using GREAT

rGREAT is an R client written by the same author of ComplexHeatmap for the web GREAT Tool


#library(devtools)
#install_github("jokergoo/rGREAT")
library(rGREAT)

# GBM.specific.H3K27ac.peaks is a GRanges object

job<- submitGreatJob(GBM.specific.H3K27ac.peaks)
tb<- getEnrichmentTables(job)
names(tb)
head(tb[[1]])
head(tb[[2]])

## job object will get updated after running getEnrichmentTables
job

Other ontologies and categories avaialble:


availableOntologies(job)
availableCategories(job)

tb<- getEnrichmentTables(job, ontology = c("GO Biological Process", "BioCyc Pathway"), category = c("GO"))
tb<- getEnrichmentTables(job, category = c("GO"))


biological.process<- tb[[1]]
biological.process[1:10,"Binom_Raw_PValue"]
head(biological.process)

mimic the enrichment bar graph:

library(ggplot2)
ggplot(biological.process[1:10,]) + 
        geom_bar(aes(x=name, y=-log10(Binom_Raw_PValue)), stat="identity", fill="blue") +
        coord_flip()

Need to re-order the name so that the y-value is ordered from high to low:


new_order<- order(-log10(biological.process$Binom_Raw_PValue), decreasing = FALSE)
biological.process$name<- factor(biological.process$name, levels = biological.process$name[new_order])

ggplot(biological.process[1:10,]) + 
        geom_bar(aes(x=name, y=-log10(Binom_Raw_PValue)), stat="identity", fill="blue") +
        coord_flip()

# or use reorder
ggplot(biological.process[1:10,]) + 
        geom_bar(aes(x=reorder(name, -Binom_Raw_PValue), y=-log10(Binom_Raw_PValue)), stat="identity", fill="blue") +
        coord_flip()

Association between genomic regions and genes can be get by plotRegionGeneAssociationGraphs(). The function will make the three plots which are same as on GREAT website and returns a GRanges object which contains the gene-region associations.

Which genes are associated with which pathway.

par(mfrow = c(1, 3))
res<- plotRegionGeneAssociationGraphs(job, ontology = c("GO Biological Process"), termID="GO:0001570")
plotRegionGeneAssociationGraphs(job)

res1<- plotRegionGeneAssociationGraphs(job, ontology = c("GO Biological Process"), termID="GO:0010574")
View(res)

crazyhottommy/rGREAT.md

Pathway analysis using GREAT