15.GO_enrichment.Rmd

---
title: "Gene Ontology Enrichment analysis"
output: 
  github_document:
bibliography: bibliography.bib
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, warning = FALSE, message = FALSE, fig.retina = 2)
library(tidyverse)
library(RColorBrewer)
require(topGO)
require(ggsci)
library(aplot)
library(ape)
library(GO.db)
library(ggtree)
library(ggstance)
library(cowplot)
```


```{r}
# If not present run 09.annotate
uniprot_gene_annot <- read_tsv("data/hpc/annotation/uniprot_gene_annot.tsv") %>% 
  mutate(go=go_ipr) %>% 
  dplyr::select(-go_ipr)

# If not present run 10.identify
candidate_regions_table <- read_rds("data/r_essentials/candidate_regions_genes_ehh.rds")

candidate_genes_table <- candidate_regions_table %>% 
  separate_rows(genes,sep=";") %>% 
  left_join(uniprot_gene_annot, by = c("genes"="geneid") ) %>% 
  unite("region_id",chr,start,end,sep = "_",remove = FALSE)
```




In a previous analysis (see [10.identify_selective_genomic_windows](10.identify_selective_genomic_windows.md) we identified a total of `r nrow(candidate_regions_table)` candidate selective sweep regions.  A summary of sweep regions and the number of genes covered by them in each population is shown below;

```{r}
candidate_genes_table %>% 
#  mutate(region_id = paste(chr,start,end,collapse = "_")) %>% 
  group_by(pop) %>% 
  dplyr::summarise(numgenes = n(), nregions = length(unique(region_id))) %>% 
  knitr::kable()
```

We then used `topGO` to identify GO terms that were enriched in the gene sets covered by selective sweeps in each population. Enrichment was initially calculated in a gene centric manner because GO terms are most properly interpreted as attached to genes not regions, however, since genomic regions and not genes are the independent units in this analysis we then recalculated enrichment statistics in a region centric manner for all GO terms found to be enriched in the initial analysis.  

The gene centric analysis considered each gene as an independent entity with the target set as all genes intersecting with selective sweeps in a given population, and the background set taken as all annotated genes for *A. digitifera*.  

To perform the region centric analysis we first assigned GO terms to regions by taking the union of all GO terms assigned to genes within a region. We then considered the target set as all candidate sweeps for a population and the background set as the complete set of all 50kb regions on which EHH statistics were calculated.

Enrichment p-values for both gene-centric and region-centric analyses were calculated using Fisher's exact test.  For the gene-centric analysis this was done using TopGO incorporating topGO's weighting system.  For the region-centric analysis it was performed using the `fisher.test` function in R.

```{r define-functions-regions}
# If not present run 10.identify
regions_genes_go <- read_rds("cache/all_regions_genes.rds") %>% 
  unite("region_id",chr,start,end) %>% 
  left_join(uniprot_gene_annot,by="geneid")

regions2go <- regions_genes_go %>% 
  group_by(region_id) %>%
  dplyr::summarise(go = paste(go,sep=";")) %>% 
  split(.$region_id) %>% 
  purrr::map(~str_split(.$go,";") %>% pluck(1) %>% str_trim() %>% unique())
```


```{r define-functions-genes}
gene2go <- uniprot_gene_annot %>% 
  split(.$geneid) %>% 
  purrr::map(~str_split(.$go,";") %>% pluck(1) %>% str_trim())

go2gene <- uniprot_gene_annot %>% 
  dplyr::select(geneid,go) %>% 
  separate_rows(go,sep=";") %>% 
  mutate(go=str_trim(go)) %>% 
  group_by(go) %>%
  split(.$go) %>% 
  purrr::map( ~.$geneid)

go_enrich_topgo <- function(targets,onto) {
    genenames <- names(gene2go)
    genelist <- factor(as.integer(genenames %in% targets))
    names(genelist) <- genenames

    GOdata <- new("topGOdata",
                  allGenes = genelist,
                  annot=annFUN.gene2GO,
                  gene2GO=gene2go,
                  ontology=onto,
                  nodeSize=10)
    
    resultFisher <- runTest(GOdata, 
                        algorithm = "weight01",
                        statistic = "fisher")

    gt <- GenTable(GOdata,classic=resultFisher,
                   WeightFisher=resultFisher,
                   orderBy = "WeightFisher",
                   topNodes=150)
    list(godata= GOdata, result=resultFisher, table=gt)
}

genes_in_enrichedGO <- function(topgo_results,target) {
  topgo_results$table %>% 
    mutate(classic=as.numeric(classic),
           WeightFisher=as.numeric(WeightFisher)) %>%
    dplyr::select(GO.ID) %>% pull %>% 
    map_df(~genesInTerm(object = topgo_results$godata,whichGO = .x)) %>%
    pivot_longer(cols = everything(),names_to = "go",values_to = "geneid") %>% 
    na.omit() %>% 
    filter(geneid %in% target)
}

go_enrich_all_ontologies <- function(genes_table,target_pop){
  target <- genes_table %>% 
    filter(pop==target_pop) %>% 
    dplyr::select(genes) %>% 
    distinct %>% 
    pull

  mf <- go_enrich_topgo(targets = target,onto = "MF")
  bp <- go_enrich_topgo(targets = target,onto = "BP")
  cc <- go_enrich_topgo(targets = target,onto = "CC")

  all_results <- list("MF"=mf,"BP"=bp,"CC"=cc)
  all_results_genes <- map_dfr(all_results,genes_in_enrichedGO,target) %>% 
    left_join(uniprot_gene_annot %>%
                dplyr::select(geneid,uniprot_id,protein)) %>%
    group_by(go) %>% 
    dplyr::summarise(uniprot_ids = paste(unique(uniprot_id),collapse = ";"), 
                geneids = paste(unique(geneid),collapse = ";"),
                protein_names = paste(unique(protein),collapse = ";"))
  
  all_onto <- rbind(bp$table %>% mutate(classic=as.numeric(classic),WeightFisher=as.numeric(WeightFisher)) %>% add_column(ontology="BP"), 
                    mf$table %>% mutate(classic=as.numeric(classic),WeightFisher=as.numeric(WeightFisher)) %>% add_column(ontology="MF"),
                    cc$table %>% mutate(classic=as.numeric(classic),WeightFisher=as.numeric(WeightFisher)) %>% add_column(ontology="CC"))
  all_onto %>% 
    left_join(all_results_genes,by=c("GO.ID"="go")) 
}
```


```{r}
if ( !file.exists("cache/topgo_results_all.rds") ){
  inshore_topgo <- go_enrich_all_ontologies(candidate_genes_table,"inshore") %>% add_column(pop="inshore")
  northoffshore_topgo <- go_enrich_all_ontologies(candidate_genes_table,"northoffshore") %>% add_column(pop="northoffshore")
  southoffshore_topgo <- go_enrich_all_ontologies(candidate_genes_table,"southoffshore") %>% add_column(pop="southoffshore")

  topgo_results_all <- rbind(inshore_topgo,northoffshore_topgo,southoffshore_topgo)
  write_rds(topgo_results_all,"cache/topgo_results_all.rds")
  write_tsv(topgo_results_all,"cache/topgo_results_all.tsv") # This table should be in Supp Info

  } else {
  topgo_results_all <- read_rds("cache/topgo_results_all.rds")
}
```


```{r}
jaccard <- function(a, b) {
    1-(length(intersect(a, b))/length(union(a,b)))
}


jaccard_go <- function(terms1){
  jaccard(terms1[[1]],terms1[[2]])
}

display_go <- topgo_results_all %>% 
  filter(Significant>0) %>% 
  pull(GO.ID) %>% 
  unique()

go2gene_signif <- topgo_results_all %>%
  filter(GO.ID %in% display_go) %>%
  dplyr::select(GO.ID,geneids) %>%
  separate_rows(geneids,sep=";") %>% 
  filter(!is.na(geneids)) %>% 
  filter(geneids!="NA") %>% 
  group_by(GO.ID) %>% 
  dplyr::summarise(geneids=paste(unique(geneids),collapse = ";")) %>% 
  ungroup() %>% 
  split(.$GO.ID) %>%
  purrr::map(~ .$geneids %>% str_split(";") %>% pluck(1))



if ( !file.exists("cache/jdm.rds")){
  cx <- cross2(go2gene_signif,go2gene_signif) #
  jdm <- cx %>% 
    map_dbl(jaccard_go) %>% 
    matrix(nrow=length(go2gene_signif))
  write_rds(jdm,"cache/jdm.rds")
} else {
  jdm <- read_rds("cache/jdm.rds")
}

rownames(jdm) <- names(go2gene_signif)
colnames(jdm) <- names(go2gene_signif)
```

```{r, fig.height=12}
n_regions_per_term <- function(term,focal_pop){
  candidate_regions_table %>%
  unite("region_id",chr,start,end,remove=FALSE) %>% 
  separate_rows(genes,sep=";") %>% 
  filter(genes %in% (go2gene_signif[term] %>% pluck(1))) %>% 
    filter(pop==focal_pop) %>% 
    pull(region_id) %>% n_distinct()
}

n_genes_per_term <- function(term,focal_pop){
  candidate_regions_table %>%
    filter(pop==focal_pop) %>% 
    separate_rows(genes,sep=";") %>% 
    filter(genes %in% (go2gene_signif[term] %>% pluck(1))) %>% 
    pull(genes) %>% n_distinct()
}

n_regions_bg_per_term <- function(term){
  regions_genes_go %>% 
    dplyr::select(region_id,geneid) %>% 
    filter(geneid %in% (go2gene_signif[term] %>% pluck(1))) %>% 
    n_distinct()
}

enrich_test <- function(n_regions,n_regions_sig,n_regions_ann,n_regions_ann_sig){
  dat <- data.frame(
  "non_signif" = c(n_regions_ann-n_regions_ann_sig, n_regions-n_regions_ann),
  "signif" = c(n_regions_ann_sig, n_regions_sig-n_regions_ann_sig),
  row.names = c("Annotated", "Non-Annotated"),
  stringsAsFactors = FALSE
  )
  colnames(dat) <- c("Non_Sig", "Sig")
  test <- fisher.test(dat)
  test$p.value
}


number_of_sig_regions <- candidate_regions_table %>% 
  group_by(pop) %>% 
  dplyr::summarise(n_sig = n()) %>% 
  pull(n_sig,name = pop)

total_bg_regions <- regions_genes_go$region_id %>% n_distinct()
```

```{r supp-table}
add_region_info <- function(data){
  data %>% 
    mutate(n_regions = map2_int(GO.ID,pop,n_regions_per_term),
         n_region_bg = map_int(GO.ID,n_regions_bg_per_term),
         n_region_expected = number_of_sig_regions[pop]*n_region_bg/total_bg_regions) %>% 
    mutate(p.region = pmap_dbl(list(total_bg_regions,number_of_sig_regions[pop],n_region_bg,n_regions),enrich_test)) 
}

supp_table6A <- topgo_results_all %>% 
  filter(classic <0.005) %>% 
  add_region_info() %>%   
  filter(n_regions>=2) %>% 
  dplyr::select(-uniprot_ids,-protein_names,-n_region_expected) %>% 
  dplyr::rename("Num Sweeps"=n_regions,"Num Background Regions"=n_region_bg,"Region p-value"=p.region)

  
supp_table6B <- supp_table6A %>% 
  separate_rows(geneids,sep = ";") %>% 
  left_join(uniprot_gene_annot,by=c("geneids"="geneid")) %>% 
  dplyr::select(GO.ID,Term,ontology,pop,geneid=geneids,uniprot_id,"Protein name"=protein)
  
write_tsv(supp_table6A,"cache/supp_table6a.tsv")
write_tsv(supp_table6B,"cache/supp_table6b.tsv")
```


```{r, fig.height=12}

bar_go <- supp_table6A %>% pull(GO.ID) %>% unique()

pg_go1 <- topgo_results_all %>% 
  filter(GO.ID %in% bar_go) %>% 
  add_region_info() 

pg_go <- pg_go1 %>% 
  pivot_wider(id_cols = c("GO.ID","ontology"), names_from = pop,values_from = p.region,values_fill = NA) %>% 
  pivot_longer(c(inshore,northoffshore,southoffshore),names_to = "pop", values_to = "p.region") %>% 
  left_join(pg_go1) %>% 
  rowwise() %>% 
  mutate(Significant = ifelse(is.na(Significant), n_genes_per_term(GO.ID,pop),Significant )) %>% 
  mutate(n_regions = ifelse(is.na(n_regions), n_regions_per_term(GO.ID,pop),n_regions)) %>% 
  mutate(n_region_bg = ifelse(is.na(n_region_bg), n_regions_bg_per_term(GO.ID), n_region_bg)) %>% 
  mutate(p.region = ifelse(is.na(p.region), enrich_test(total_bg_regions, number_of_sig_regions[pop],n_region_bg, n_regions), p.region)) 
```




We found a total of `r length(bar_go)` GO terms that were enriched (gene and region centric p<0.005; at least 2 regions) in at least one of the three of the locations.  These are summarised in Figure 1.  


```{r, fig.height=12}
pop_names <- c("inshore"="Inshore","northoffshore"="North Offshore","southoffshore"="South Offshore")

library(colorspace)
hclcols <- sequential_hcl(3,palette = "ag_Sunset")

# Text labels have separate data so we can show them when bars are absent
text_label_data <- pg_go %>% 
  mutate(text_label = paste(Significant,"/",n_regions,sep = "")) %>% 
  mutate(pop_names=pop_names[pop]) %>% 
  mutate(classic = case_when(
    classic<0.1 ~ classic,
    classic>=0.1 ~ 0.8
  ))

barplot <- pg_go %>% 
  filter(classic<0.1) %>% 
  mutate(pop_names=pop_names[pop]) %>% 
  mutate(text_label = paste(Significant,"/",n_regions,sep = "")) %>% 
  mutate(is_enriched = ifelse(p.region<0.005, 1, 0.3)) %>% 
  ggplot(aes(x=-log10(classic),y=GO.ID)) + 
  geom_col(aes(alpha=is_enriched,fill=ontology,color=ontology),size=0.1) +
  scale_color_manual(values = hclcols,labels = c("Biological process", "Cellular component", "Molecular function")) +
  scale_fill_manual(values = hclcols,labels = c("Biological process", "Cellular component", "Molecular function")) +  
  geom_text(data=text_label_data,aes(label=text_label), hjust=-0.1,size=2.5) +
  theme_bw() + 
  guides(alpha="none") +
  labs(x=expression(-Log[10](P)),y="",fill="") +
  theme(panel.grid = element_blank(),
        legend.position = "bottom",
        legend.text = element_text(),
        axis.text.y = element_blank(),
        axis.text.x = element_text(size=9),
        axis.title.x = element_text(size=11),
        strip.text = element_text(size=9),
        strip.background = element_blank()) + 
  ylab("") + 
  facet_wrap(~pop_names) +
  xlim(0,7)

leg_b <- get_legend(barplot + guides(alpha="none",color="none"))

go_terms = AnnotationDbi::select(GO.db, keys(GO.db, "GOID"), c("TERM", "ONTOLOGY")) %>% 
  filter(GOID %in% display_go)



tree_data <- pg_go %>% 
  pivot_wider(id_cols = c("GO.ID","ontology"), names_from = pop, values_from = classic) %>% 
  dplyr::rename(inshore_p = inshore, northoffshore_p = northoffshore, southoffshore_p = southoffshore) %>% 
  left_join(pg_go %>% pivot_wider(id_cols = c("GO.ID"), names_from = pop, values_from = Significant)) %>% 
  as.data.frame() %>% 
  left_join(go_terms,by=c("GO.ID"="GOID")) %>% 
  mutate(term_label = paste(GO.ID," ",TERM,sep=""))

jdm_included <- jdm[tree_data$GO.ID,tree_data$GO.ID]
hc <- hclust(as.dist(jdm_included))
hc_order <- hc$order
names(hc_order) <- rownames(jdm_included)[hc$order]

# Edit term labels to make them shorter
tree_data_trim <- tree_data %>% 
  mutate(term_label = case_when(
    term_label=="GO:0000981 DNA-binding transcription factor activity, RNA polymerase II-specific" ~ "GO:0000981 DNA-binding transcription factor activity",
    term_label=="GO:0007186 G protein-coupled receptor signaling pathway" ~ "GO:0007186 G protein-coupled receptor signaling",
#    term_label=="GO:0016879 ligase activity, forming carbon-nitrogen bonds" ~ "GO:0016879 ligase activity, forming C-N bonds",
    TRUE ~ term_label
  ))

tree <- ggtree(as.phylo(hc)) %<+% tree_data_trim + 
geom_tiplab(aes(label = term_label),align = TRUE, size=2.8) + 
  xlim(0,10) 

bpt <- (barplot + theme(legend.position = "none")) %>% insert_left(tree, width = 1.2)
#ggsave(bpt,filename = "figures/go_enrichment_wide.png",height = 2,width = 6.75)
ggsave(bpt,filename = "figures/go_enrichment_wide.pdf",height = 5.5,width = 16.9, units = "cm")
leg_b <- ggdraw(leg_b)
ggsave(leg_b,filename = "figures/go_legend.pdf",width = 16.9,height = 1.25,units = "cm")


bpt
```

**Figure 1: Enriched GO terms for genes intersecting with candidate loci under selection**. Length of bars indicates significance of gene-centric enrichment (longer is more significant) and colour indicates the ontology. Filled bars are also significant (p<0.005) in a region-centric test while unfilled bars are not. Numerical labels are given as ng/nr where ng indicates the number of genes contained within sweep regions and nr indicates the number of distinct regions where these genes are found. Dendrogram depicts relationships between GO terms based on numbers of shared annotated genes.