library(topGO)

# Read in GO mappings to gene ids
# <gene>\t<GOID>,<GOID>,...
vng2GO <- readMappings(file='vng2GO.map')

# Load up cluster or gene set identifiers
geneSets <- ## Depends on your analysis, list of vectors of gene ids

################
# Run analysis #
################
#  Make Gene list
geneNames <- # Need to detemine what this is
geneList <- factor(as.integer(geneNames %in% tmp1))
names(geneList) <- geneNames
# Make Biological Process GOData object
tmp1 <- c('VNG####','VNG####') # put in a couple of gene ids to build structure for automation
GOdata.BP <- new("topGOdata", ontology='BP', allGenes = geneList, annot = annFUN.gene2GO, gene2GO = vng2GO)
m1.BP <- matrix(nrow = length(GOdata.BP@graph@nodes), ncol = length(geneSets), dimnames = list(GOdata.BP@graph@nodes, 1:length(geneSets)))
m2.BP <- matrix(nrow = length(geneSets), ncol=2, dimnames = list(1:length(geneSets), c('Top10.Terms.BP', 'BH.sig.GO.Ids.BP')))

# Make Molecular Function GOData object
#GOdata.MF <- new("topGOdata", ontology='MF', allGenes = geneList, annot = annFUN.gene2GO, gene2GO = vng2GO)
#m1.MF <- matrix(nrow = length(GOdata.MF@graph@nodes), ncol = length(geneSets), dimnames = list(GOdata.MF@graph@nodes, 1:length(geneSets)))
#m2.MF <- matrix(nrow = length(geneSets), ncol=2, dimnames = list(1:length(geneSets), c('Top10.Terms.MF', 'BH.sig.GO.Ids.MF')))

# Make Cellular Component GOData object
#GOdata.CC <- new("topGOdata", ontology='CC', allGenes = geneList, annot = annFUN.gene2GO, gene2GO = vng2GO)
#m1.CC <- matrix(nrow = length(GOdata.CC@graph@nodes), ncol = length(geneSets), dimnames = list(GOdata.CC@graph@nodes, 1:length(geneSets)))
#m2.CC <- matrix(nrow = length(geneSets), ncol=2, dimnames = list(1:length(geneSets), c('Top10.Terms.CC','BH.sig.GO.Ids.CC')))

for( geneSet in (1:length(geneSets)) ) {
    # Expand gene list and change factor in GOdata
    genes <- geneSets[[geneSet]]
    GOdata.BP@allScores <- factor(as.integer(geneNames %in% genes))
    #GOdata.MF@allScores <- factor(as.integer(geneNames %in% genes))
    #GOdata.CC@allScores <- factor(as.integer(geneNames %in% genes))
    # Biological process
    r1.BP = runTest(GOdata.BP, algorithm = 'classic', statistic = 'fisher')
    m1.BP[,biclust] = r1.BP@score
    m2.BP[biclust,1] = paste(GenTable(GOdata.BP, r1.BP)[,2],collapse=', ')
    m2.BP[biclust,2] = paste(names(which(p.adjust(r1.BP@score,method='BH')<=0.05)),collapse=', ')
    #r1.MF = runTest(GOdata.MF, algorithm = 'classic', statistic = 'fisher')
    #m1.MF[,biclust] = r1.MF@score
    #m2.MF[biclust,1] = paste(GenTable(GOdata.MF, r1.MF)[,2],collapse=', ')
    #m2.MF[biclust,2] = paste(names(which(p.adjust(r1.MF@score,method='BH')<=0.05)),collapse=', ')
    #r1.CC = runTest(GOdata.CC, algorithm = 'classic', statistic = 'fisher')
    #m1.CC[,biclust] = r1.CC@score
    #m2.CC[biclust,1] = paste(GenTable(GOdata.CC, r1.CC)[,2],collapse=', ')
    #m2.CC[biclust,2] = paste(names(which(p.adjust(r1.CC@score,method='BH')<=0.05)),collapse=', ')
}

# Write out the results, can add the others by doing cbind to add the other matrices as columns
write.csv(m2.BP,'biclusterEnrichment_GOBP.csv')