#********************************************************************************************************
# Build a customized object of type topGOdata for Paramecium macronuclear genes v2
#********************************************************************************************************
#source("https://bioconductor.org/biocLite.R")
#biocLite('topGO')
#biocLite("GO.db")
#biocLite('ALL')
#biocLite('hgu95av2.db')
library(topGO)
# Path to file_all Paramecium genes mapped to go terms (BP or MF)
Path_to_genes_GO_terms <- " path to Coding-Genes_BP_v2.0.txt "
# Path to file_IES-containing genes mapped to go terms (BP or MF)
Path_to_IES_genes_GO_terms <- " path to IES-Genes_BP_v2.0.txt "

# Build geneID2GO object from the desired gene universe Gene_ID:GO_MF/GO_BP (for all 40461 macronuclear coding gene models from v2.0)
geneID2GO <- readMappings(Path_to_genes_GO_terms, sep = "\t")
str(head(geneID2GO))
geneNames <- names(geneID2GO[-1])
head(geneNames)

## Create a set of lists of interesting genes (e.g.IES-containing genes) 
# All Paramecium 19118 IES-containing genes
myInterestingGenes <- read.table(Path_to_IES_genes_GO_terms, sep = "\t")
colnames(myInterestingGenes) <- "Gene_ID"
geneList <- factor(as.integer(geneNames%in%myInterestingGenes$Gene_ID))
names(geneList) <- geneNames
str(geneList)

GOdataMF <- new("topGOdata",
                description = "Functional Enrichment Analysis of IES-containing genes", ontology = "BP", # Change to BP/MF
                allGenes = geneList, 
                nodeSize = 5,
                annot = annFUN.gene2GO, gene2GO = geneID2GO)

GOdataMF

# Enrichment Tests_Fisher for count data
# Test the over-representation of GO terms within the group of interesting genes
resultFisher <- runTest(GOdataMF, statistic = "fisher", algorithm = "weight01")
resultFisher

# Analysis of Results
# GeneTable
allGO = usedGO(object = GOdataMF)
allRes <- GenTable(GOdataMF, classicFisher = resultFisher,
                   orderBy = "classicFisher", ranksOf = "classicFisher", topNodes =length(allGO))

# Adjustment for multiple comparison
allRes$padjust <- p.adjust(allRes$classicFisher, method = "hochberg")
allRes_adj <-allRes[allRes$padjust < 0.05,]
allRes_adj

# Write results to table
#write.csv(allRes_adj, "/Users/cataniastudent/Desktop/GO-IESall_Table_BP.txt", row.names = F, quote = F)

