Methods_S4.Rmd

```{r}
#METHODS S4 GEne ontology analysis for molecular mechanisms
```

```{r}

library(grex);
library(biomaRt);
library(topGO);
library(dplyr);
library(data.table);
library(ggplot2)

```

```{r}
annot = fread("02_GSandMM.csv",
              select = c("geneSymbol","ENSG","moduleColor"));
annot$ensembl = substr(annot$ENSG,1,15);
ensembl_id = annot$ensembl;
all_id = grex(ensembl_id);

annotComplete = merge(annot, all_id,
                     by.x="ensembl",
                     by.y = "ensembl_id");

write.csv(annotComplete, "03_moduleGeneAnnotated.csv")
```

```{r}
genes_bg = annotComplete$entrez_id;
length(genes_bg)
```

```{r}
tot_background = length(genes_bg);
db = useMart("ENSEMBL_MART_ENSEMBL",
             dataset= "hsapiens_gene_ensembl",
             host = "https://www.ensembl.org");
go_ids = getBM(attributes = c("go_id",
                              "entrezgene_id",
                              "namespace_1003"),
               filters = "entrezgene_id",
               values = genes_bg,
               mart=db)
```


```{r}
modInt = as.factor(annotComplete$moduleColor);
annotSplit = split(annotComplete, modInt);
candidate_list = annotSplit$yellow$entrez_id;

length(candidate_list);

```

```{r}
tot_candidate = length(candidate_list);
keep = candidate_list %in% go_ids[,2];
keep = which(keep ==TRUE);
candidate_list = candidate_list[keep];
geneList = factor(as.integer(genes_bg %in% candidate_list));
names(geneList) = genes_bg
```

```{r}
gene2GO = unstack(go_ids[,c(1,2)]);

GOdata = new("topGOdata",
             ontology = c("BP"),
             allGenes = geneList,
             annot = annFUN.gene2GO,
             gene2GO = gene2GO);

allGO = usedGO(GOdata)
```

```{r}
sigTest = runTest(GOdata,
                  algorithm = "elim",
                  statistic = "fisher")
```

```{r}
all_res = GenTable(GOdata,
                   weightFisher = sigTest,
                   orderBy = "weightFisher",
                   topNodes = length(allGO))
```

```{r}
all_res$OR = log2((all_res$Significant/tot_candidate)/(all_res$Annotated/tot_background));
write.csv(all_res,"03_moduleGOAnnotated.csv")
```

```{r}
GO_bar = all_res %>%
  filter(weightFisher<0.05)%>%
  filter(Annotated>=30);
GO_bar = dplyr::select(GO_bar,Term,weightFisher,OR);
sapply(GO_bar,class);

```

```{r}
GO_bar = transform(GO_bar,weightFisher = as.numeric(weightFisher));
sapply(GO_bar,class);
```

```{r}
GO_bar = GO_bar %>%
  dplyr::arrange(weightFisher)%>%
  dplyr::slice(1:10);
CapStr = function(y){
  c = strsplit(y," ")[[1]]
  paste(toupper(substring(c,1,1)),substring(c,2),
        sep="",collapse=" ")
}

GO_bar$Term = sapply(GO_bar$Term,CapStr);

p =ggplot(GO_bar,
          aes(x=reorder(Term,OR),y=OR,fill=weightFisher)) +
  geom_bar(stat = "identity", width =0.9, color = "black") +
  coord_flip() +
  scale_fill_gradient(low="#feff2b",high = "#fe0100") +
  ylim(0,3)+
  labs(title = ~underline("Enriched GO Biological Processes"),
       x = NULL,
       y = "Odds Ratio",
       fill = "p.value") +
  theme_bw() +
  theme(plot.title.position = "plot") +
  theme(plot.title = element_text(size=12))+
  theme(axis.title.x=element_text(size=10,face="bold"),
        axis.title.y = element_text(size=10,face="bold")) +
  theme(axis.text.x = element_text(size=10, face = "bold"),
        axis.text.y = element_text(size=10,face = "bold")) +
  theme(legend.position = "right");
p
```