absREL.GOenrichment.Rmd

---
title: "Inquiline Genomics: GO enrichment analysis of absREL results"
output:
  html_document:
    toc: true
    toc_float: true
    toc_depth: 2
    collapsed: true
    number_sections: true
    df_print: paged
  pdf_document: default
---
<STYLE TYPE="text/css">
<!--
  td{
    font-family: Arial; 
    font-size: 8pt;
    padding:0px;
    cellpadding="0";
    cellspacing="0"
  }
  th {
    font-family: Arial; 
    font-size: 8pt;
    height: 20px;
    font-weight: bold;
    text-align: right;
    background-color: #ccccff;
  }
  table { 
    border-spacing: 0px;
    border-collapse: collapse;
  }
--->
</STYLE>





# R steps
```{r}
source("https://bioconductor.org/biocLite.R")
BiocManager::install("GO.db")
biocLite("GO.db")
library("topGO")
library(dplyr)
library(plyr)
library(tidyr)
```

### Read large input table
```{r}
all<-read.csv("/Users/lukas/sciebo/inquilineGenomics18/SOS/data/absREL.tsv",sep="\t")
```

### create topGO element for Cluster GO annotations
```{r}

scoID2go<-list()
  for(i in 1:length(all$OGid)){
  scoID2go[[i]]<-unlist(strsplit(as.character(all$GO[i]),", "))
  names(scoID2go)[i]<-as.character(all$OGid[i])
  }

```

### Prepare Objects for GO enrichment analysis
```{r}

scoNames.all<-names(scoID2go)
#restrict gene space to all SCOs orthologs
scoNames<-scoNames.all
scoNames <- factor(scoNames)
names(scoNames) <- scoNames

```


### Select Clusters of Interest from dataframe "all".
```{r}
# create a 0/1 matrix (interesting/not interesting)
######################################

## all parasites
#SCOOfInterest<-unlist(subset(all,p.adjust(P1pv,method = "fdr")<0.05 | p.adjust(AINSpv,method = "fdr")<0.05 | p.adjust(PARGpv,method = "fdr")<0.05 | p.adjust(ACHApv,method = "fdr")<0.05,OGid))

# Origin of parasitism
#SCOOfInterest<-unlist(subset(all,PARGfdr<0.05 | AINSfdr<0.05 | ACHAfdr<0.05 |  P1fdr<0.05,OGid))
SCOOfInterest<-unlist(subset(all,p.adjust(all$AECHpv,method = "fdr")<0.1 | p.adjust(all$AHEYpv,method = "fdr")<0.1,OGid))
#SCOOfInterest<-unlist(subset(all,AINSfdr<0.1 |  P1fdr<0.1,OGid))
#SCOOfInterest<-unlist(subset(all,PARGfdr<0.05,OGid))

#SCOOfInterest<-unlist(subset(all,p.adjust(PARGpv,method = "fdr")<0.05 ,OGid))
#SCOOfInterest<-unlist(subset(all,p.adjust(ACHApv,method = "fdr")<0.05 ,OGid))


scoList <- factor(as.integer(scoNames %in% SCOOfInterest))

table(scoList)
names(scoList) <- scoNames

posSelSet<-subset(all,OGid %in% SCOOfInterest)

```


```{r}
# create object for Biological Process, Molecular Function and Cellular component:
##################################################################################

GOdata_MF <- new("topGOdata", ontology = "MF", allGenes = scoList, annot = annFUN.gene2GO, gene2GO = scoID2go)
GOdata_BP <- new("topGOdata", ontology = "BP", allGenes = scoList, annot = annFUN.gene2GO, gene2GO = scoID2go)
GOdata_CC <- new("topGOdata", ontology = "CC", allGenes = scoList, annot = annFUN.gene2GO, gene2GO = scoID2go)

# Molecular Function
####################
resultFis_MF 	<- runTest(GOdata_MF, algorithm = "parentchild", statistic = "fisher")
table_MF <- GenTable(GOdata_MF, parentChild = resultFis_MF, orderBy = "parentChild", ranksOf = "parentChild", topNodes = 200)

# Biological Process
####################
resultFis_BP <- runTest(GOdata_BP, algorithm = "parentchild", statistic = "fisher")
table_BP <- GenTable(GOdata_BP, parentChild = resultFis_BP, orderBy = "parentChild", ranksOf = "parentChild", topNodes = 200)


# Cellular Component
####################
resultFis_CC <- runTest(GOdata_CC, algorithm = "parentchild", statistic = "fisher")
table_CC <- GenTable(GOdata_CC, parentChild = resultFis_CC, orderBy = "parentChild", ranksOf = "parentChild", topNodes = 200)


# filter out significant GO terms:
###################################

MF <- subset(table_MF, parentChild < 0.1)
BP <- subset(table_BP, parentChild < 0.1)
CC <- subset(table_CC, parentChild < 0.1)

MF$ontology<-"MF"
BP$ontology<-"BP"
CC$ontology<-"CC"

MF
BP
CC
```
```{r}

subset(all,OGid %in% genesInTerm(GOdata_BP, "GO:0006928")[[1]] & (AINSfdr<0.05 | PARGfdr<0.05|P1fdr<0.05|ACHAfdr<0.05))
subset(all,OGid %in% genesInTerm(GOdata_BP, "GO:0006030")[[1]] & (AINSfdr<0.05 | PARGfdr<0.05|P1fdr<0.05|ACHAfdr<0.05))
subset(all,OGid %in% genesInTerm(GOdata_BP, "GO:0023051")[[1]] & (AINSfdr<0.05 | PARGfdr<0.05|P1fdr<0.05|ACHAfdr<0.05))

subset(all,OGid %in% genesInTerm(GOdata_BP, "GO:0000272")[[1]] & (AINSfdr<0.05 | PARGfdr<0.05|P1fdr<0.05|ACHAfdr<0.05))

subset(all,OGid %in% genesInTerm(GOdata_MF, "GO:0008289")[[1]] & (AINSfdr<0.05 | PARGfdr<0.05|P1fdr<0.05|ACHAfdr<0.05))



```

### save results to file
```{r}
# Origins of social parasitism
res<-rbind(MF,BP) %>% rbind(.,CC)
outfile<-"/Users/lukas/sciebo/inquilineGenomics18/GOenrichment/results/GOenrichment.absREL.parasiteOrigins.tsv"
#outfile<-"/Users/lukas/sciebo/inquilineGenomics18/GOenrichment/results/GOenrichment.absREL.parasiteBranches.tsv"
#write.table(res,outfile,quote=F,sep="\t",row.names = F)
cat("# GO enrichment of genes under positive selection at either social parasite origin (AINS or P1) (at FDR<0.1). \n",file=outfile)
#cat("# GO enrichment of genes under positive selection at either social parasite branch (ACHA, PARG, AINS or P1) (at FDR<0.1). \n",file=outfile)
#write.table(res, outfile,append=TRUE, quote=F,sep="\t",row.names=F)


```