---
title: "EnrichResults in network"
author: "Carine"
date: "2025"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(tidyverse)
```


# Load data

```{r}
total_DE <- read_tsv("Enrichment/GO_clusterProfiler_DE_total.tsv") %>% filter(ONTOLOGY == "BP")
total_WGCNA <- read_tsv("Enrichment/GO_clusterProfiler_WGCNA_total.tsv")  %>% filter(ONTOLOGY == "BP")
dim(total_WGCNA)
recent_DE <- read_tsv("Enrichment/GO_clusterProfiler_DE_recent.tsv")  %>% filter(ONTOLOGY == "BP")
recent_WGCNA <- read_tsv("Enrichment/GO_clusterProfiler_WGCNA_recent.tsv")  %>% filter(ONTOLOGY == "BP")

ancient_DE <- read_tsv("Enrichment/GO_clusterProfiler_DE_ancient.tsv") 
dim(ancient_DE)
ancient_WGCNA <- read_tsv("Enrichment/GO_clusterProfiler_WGCNA_ancient.tsv")  
dim(ancient_WGCNA)

murinae_DE <- read_tsv("Enrichment/GO_clusterProfiler_DE_murinae.tsv")  
dim(murinae_DE)
murinae_WGCNA <- read_tsv("Enrichment/GO_clusterProfiler_WGCNA_murinae.tsv")  %>% filter(ONTOLOGY == "BP")
dim(murinae_WGCNA)
```

## create the function

```{r}
library(igraph)
library(tidygraph)
library(ggraph)
#library(ggnetwork)
library(ggforce)
library(ggalt) #geom_encircle

library(ggthemes)


plot_enrichResult_table_in_network <- function(
    enrich_res_table_tmp ,
    output_path=NULL,
    width_output=6,
    height_output=6,
    title = NULL,
    resolution_clustering=1,
    threshold_quantil_edges = 0,
    threshold_percent_edges = 0,
    force_attraction_cluster=2,
    force_repulsion_name=10,
    seed_tmp=123,
    plot_edge_number=F,
    plot_GO_name = T,
    increase_plot_area=1.1,
    cluster_name_df=NULL) {
  
  if(!nrow(enrich_res_table_tmp) >0) {
    return(NULL)
  }
  
  ## create a simplified table with
  ###the list of genes per GO (geneID_l)
  ###the geneRation in numeric GeneRatio_num
  ###and a column with the GO name split in small lines (Description2)
  
  df_go_res_simplified <- enrich_res_table_tmp  %>%
    rowwise() %>%
    mutate(geneID_l = str_split(geneID, "/"),
           GeneRatio_1= ((unlist(str_split(GeneRatio,"/"))[1])),
           GeneRatio_2= ((unlist(str_split(GeneRatio,"/"))[2])),
           GeneRatio_num = as.numeric(GeneRatio_1)/as.numeric(GeneRatio_2),
           Description2 = paste0(strwrap(Description, 25), collapse="\n")
    ) %>%
    dplyr::select(Description,Description2,ID, Count,qvalue,ONTOLOGY,geneID_l,GeneRatio_num,GeneRatio,p.adjust)
  
  
  ## compute the number of shared gene between all GO
  
  df_go_res_simplified_all_combi <- as.tibble(t(combn(df_go_res_simplified$Description,2))) %>%
    dplyr::rename(gs1 = "V1", gs2 = "V2") %>%
    left_join(df_go_res_simplified%>%select(Description,geneID_l), by = c("gs1"="Description")) %>%
    left_join(df_go_res_simplified%>%select(Description,geneID_l), by = c("gs2"="Description"), suffix = c(".gs1",".gs2")) 
  
  df_go_res_simplified_all_combi <- df_go_res_simplified_all_combi %>%
    rowwise()%>%
    mutate(l_union = length(union( geneID_l.gs1,geneID_l.gs2)),
           n_common = length(intersect( geneID_l.gs1,geneID_l.gs2)),
           p_shared = n_common/l_union )  %>%
    arrange(-p_shared)
  
  ## get only interesting column
  GO_ovlap <- df_go_res_simplified_all_combi %>% ungroup() %>% dplyr::select(gs1,gs2,p_shared,n_common)
  
  #create an overlap network
  GO_ovnet <-  igraph::graph_from_data_frame(GO_ovlap,directed = F, vertices = df_go_res_simplified)
  
  
  # create weighted edges according the percentage of shared genes
  E(GO_ovnet)$weight <- edge_attr(GO_ovnet)$p_shared*10
  GO_ovnet <- delete_edges(GO_ovnet,E(GO_ovnet)[weight==0])
  
  
  set.seed(seed_tmp)
  
  # Define clusters of GO based on weighted edges # A larger edge weight means a stronger connection for this function
  cl_louvain <-cluster_louvain(GO_ovnet, resolution = resolution_clustering)
  grp_louvain <- as.factor(cl_louvain$membership)
  nb_grp_louvain <- length(levels(grp_louvain))
  df_louvain <- tibble(title =title,
                       Description = names(membership(cl_louvain)),
                       cluster = cl_louvain$membership)
  
  vertex_attr(GO_ovnet)$cluster_louvain = factor(cl_louvain$membership)
  
  
  # give other weights to build the layout (where go are in the graph) for a better separating of the cluster on the graph
  weights2 <- ifelse(crossing(cl_louvain, GO_ovnet), 1, force_attraction_cluster)
  weights3 <- weights2*E(GO_ovnet)$p_shared*5
  weights3 <- ifelse(weights3 > force_attraction_cluster, weights3, weights3)
  layoutx <- layout_with_fr(GO_ovnet, weights=weights3)
  
  # keep only a subset of the edges
  GO_ovnet_cut = GO_ovnet
  
  if(threshold_quantil_edges>0) {
    cut.off <- quantile(E(GO_ovnet_cut)$weight , threshold_quantil_edges)
    GO_ovnet_cut <- delete_edges(GO_ovnet_cut,E(GO_ovnet_cut)[weight<=cut.off])
  }
    # keep only a subset of the edges
  if(threshold_percent_edges>=0 & threshold_percent_edges <= 1) {
    GO_ovnet_cut <- delete_edges(GO_ovnet_cut,E(GO_ovnet_cut)[p_shared<=threshold_percent_edges])
  }
  
  
  
  ggobj <- GO_ovnet_cut%>% 
    ggraph(layout=layoutx) + theme_graph(base_family="", base_size = 12) 
  
  coord_data <- ggobj$data
  
  
  ## get cluster position to numerotate GO
  
   cluster_coord <- coord_data %>%
      group_by(cluster_louvain) %>%
      summarise(
        x_clust = mean(x),
        y_clust = mean(y)) %>%
     arrange(-y_clust,x_clust) %>%
    mutate(layout_cluster = 1:nb_grp_louvain)
     
      
  df_louvain_ok <- coord_data %>% left_join(df_louvain, by = c("name"="Description")) %>% 
    left_join(cluster_coord) %>%
    dplyr::select(title,x,y,name,cluster_louvain, cluster,.ggraph.orig_index,layout_cluster,ID,Count,GeneRatio,p.adjust) %>%
    arrange(layout_cluster,-y, x) %>%
    mutate(plot_ID = 1:length(grp_louvain)) %>%
    arrange(.ggraph.orig_index)
  
  
  ggobj <-ggobj +
    # Create circle around points within group and label +
    geom_encircle(aes(
      x,
      y,
      color = cluster_louvain,
      fill = cluster_louvain
    ),
    s_shape=0.2, expand=0.08, spread=0.01,
    na.rm = T,
    alpha = 0.3,
    show.legend = FALSE)
  
  
  if(plot_edge_number) {
    ggobj <- ggobj +  geom_edge_link(aes(edge_width=p_shared * 100, alpha= p_shared * 100, label = n_common),
                                     color="grey50",
                                     angle_calc = 'along',
                                     label_size = 2,
                                     label_colour = "darkblue")
  } else {
    ggobj <- ggobj +  geom_edge_link(aes(edge_width=p_shared * 100, alpha= p_shared * 100),
                                     color="grey50",
                                     angle_calc = 'along',
                                     label_size = 2,
                                     label_colour = "darkblue")
  }
  
  
  if(plot_GO_name) {
    ggobj <-  ggobj + geom_node_text(aes(label=vertex_attr(GO_ovnet)$Description2, color = cluster_louvain),
                                     repel = T, size = 3, box.padding = 1, force = force_repulsion_name,
                                     max.overlaps = 40, max.time = 4, max.iter = Inf,
                                     segment.alpha = 0.5, show.legend = F)
  } 
  
  ggobj <-  ggobj + 
    geom_node_point(aes(shape=ONTOLOGY, size = GeneRatio_num, fill=cluster_louvain)) +
    geom_node_text(data = df_louvain_ok, aes(label=plot_ID, x=x, y= y), color="black", size = 4,repel = F) 
  
  
  ggobj <-  ggobj +
    scale_shape_manual(values = c("MF"=23, "CC"=22, "BP" = 21)) + 
    scale_size_continuous(name="GeneRatio", range = c(6,9)) +
    scale_edge_width_continuous(name="% DEGs shared", breaks = seq(max(threshold_percent_edges*100,0),100, length=3),limits = c(0,100), range = c(1,3)) +
    scale_edge_alpha_continuous(name="% DEGs shared", breaks = seq(max(threshold_percent_edges*100,0),100, length=3), limits = c(0,100), range = c(0.2,0.8)) +
    scale_fill_ptol() +
    scale_color_ptol() +
    guides(fill = "none", shape = guide_legend(override.aes = list(size = 5))) + theme(legend.position = "top",  legend.box="vertical") 
  
  
  
  xlim_tmp <- ggplot_build(ggobj)$layout$panel_scales_x[[1]]$range$range 
  ylim_tmp <- ggplot_build(ggobj)$layout$panel_scales_y[[1]]$range$range 
  
  if (increase_plot_area) {
    
    xlim_tmp_new <- c(median(xlim_tmp)-abs((median(xlim_tmp) - xlim_tmp[1] )*increase_plot_area), median(xlim_tmp) +abs(xlim_tmp[2]-median(xlim_tmp))*increase_plot_area)
    ylim_tmp_new <- c(median(ylim_tmp)-abs((median(ylim_tmp) - ylim_tmp[1] )*increase_plot_area), median(ylim_tmp) +abs(ylim_tmp[2]-median(ylim_tmp))*increase_plot_area)
    
    ylim_tmp <- xlim_tmp_new
    ylim_tmp <- ylim_tmp_new
    ggobj <-  ggobj +  xlim(xlim_tmp_new) +  ylim(ylim_tmp_new)
    
    
  }
  
  if(is.character(title)) {
    ggobj <- ggobj + ggtitle(title)
  }
  #ggsave(paste(save_dir_tmp, "test_representation_cluster_fig2G_with_name.pdf", sep = ""), width =15, height = 15)
  ggsave( plot=ggobj, filename = output_path , width =width_output, height = height_output)
  
  
  df_louvain_ok %>% arrange(layout_cluster, plot_ID) %>%
    dplyr::select(title,plot_ID,cluster_louvain,ID,name,GeneRatio,p.adjust) %>%
    dplyr::rename(Description="name",
                  group_ID = "cluster_louvain") %>%
    write_tsv(str_replace(output_path,".pdf$",".tsv"))
  

  
  ############ add cluster name using cluster name
  
  ggobj2 = NULL
  if( ! is.null(cluster_name_df)) {
    
    ggobj2 <- ggobj
    ##remove go name
    if(plot_GO_name) {
      ggobj2$layers[[3]] <-NULL
    }
    
    cluster_coord <- cluster_coord  %>%  rowwise()  %>%
      mutate(cluster_louvain_num=as.numeric(cluster_louvain)) %>%
      left_join(cluster_name_df, by = c("cluster_louvain_num" = "group_ID")) %>%
      mutate(cluster_name_swap = paste0(strwrap(cluster_name, 20), collapse="\n")) %>%
      ungroup()
    
    
    xlim_tmp_new2 <- c(median(xlim_tmp)-abs((median(xlim_tmp) - xlim_tmp[1] )*2), median(xlim_tmp) +abs(xlim_tmp[2]-median(xlim_tmp))*2)
    ylim_tmp_new2 <- c(median(ylim_tmp)-abs((median(ylim_tmp) - ylim_tmp[1] )*2), median(ylim_tmp) +abs(ylim_tmp[2]-median(ylim_tmp))*2)

    
    
    ggobj2 <- ggobj2 +
      ggrepel::geom_label_repel(
        data = cluster_coord %>% filter(x_clust >= median(xlim_tmp)),
        aes(
          x = x_clust,
          y = y_clust,
          label = cluster_name_swap,
          color = cluster_louvain
        ),
        size = 5,
        segment.size = 0.2,
        show.legend = F,
        # arrow = arrow(
        #   length = unit(0.03, "npc"),
        #   type = "closed",
        #   ends = "first"
        # ),
        point.padding = NA,
        box.padding = 0.5,
        xlim  = c(xlim_tmp[2]-(xlim_tmp[2] - median(xlim_tmp))/1.5, NA)
      ) +
      ggrepel::geom_label_repel(
        data = cluster_coord %>% filter(x_clust < mean(xlim_tmp)),
        aes(
          x = x_clust,
          y = y_clust,
          label = cluster_name_swap,
          color = cluster_louvain
        ),
        size = 5,
        segment.size      = 0.2,
        show.legend = F,
        # arrow = arrow(
        #   length = unit(0.03, "npc"),
        #   type = "closed",
        #   ends = "first"
        # ),
        point.padding = NA,
        box.padding = 0.5,
        xlim  = c(NA, xlim_tmp[1]+(median(xlim_tmp)-xlim_tmp[2])/1.5)
      ) +
      xlim(xlim_tmp_new2)  +
      ylim(ylim_tmp_new2)
    
    ggsave( plot=ggobj2, filename = str_replace(output_path, ".pdf$", "_with_cluster_name.pdf") , width =width_output, height = height_output)
    
    df_louvain_ok <- df_louvain_ok %>% left_join(cluster_name_df, by=c("cluster"="group_ID")) 
    
    df_louvain_ok %>% arrange(layout_cluster, plot_ID) %>%
      dplyr::select(plot_ID,cluster_louvain,cluster_name,ID,name,GeneRatio,p.adjust) %>%
      dplyr::rename(Description="name",
                    group_ID = "cluster_louvain",
                    group_Description = "cluster_name") %>%
      write_tsv(str_replace(output_path, ".pdf$", "_with_cluster_name.tsv"))
     
  }
  
  
  
  
  
  return(list(ggobj,ggobj2) )
}

```


## Use function 

```{r}
save_dir_tmp <- "output_plot_enrichResult_in_network_final"
dir.create(save_dir_tmp)

```


### Ancient

```{r}
p1_ancient_DE <- plot_enrichResult_table_in_network(
  enrich_res_table_tmp = ancient_DE,
  output_path=paste(save_dir_tmp, "/test_representation_cluster_ancient_DE_with_GO_name.pdf", sep = ""),
  width_output=8,
  height_output=8,
  title = "ancient_DE",
  resolution_clustering=1,
  threshold_percent_edges = 0.1,
  force_attraction_cluster=2,
  force_repulsion_name=40,
  seed_tmp=123,
  plot_edge_number=F,
  plot_GO_name = T,
  increase_plot_area=1.1)
p1_ancient_DE
```

```{r, fig.height=6,fig.width=6}
p1_ancient_WGCNA <- plot_enrichResult_table_in_network(
  enrich_res_table_tmp = ancient_WGCNA,
  output_path=paste(save_dir_tmp, "/test_representation_cluster_ancient_WGCNA_with_GO_name.pdf", sep = ""),
  title = "ancient_WGCNA",
  width_output=6,
  height_output=6,
  resolution_clustering=1,
  threshold_percent_edges = 0.1,
  force_attraction_cluster=6,
  force_repulsion_name=50,
  seed_tmp=123,
  plot_edge_number=F,
  plot_GO_name = T,
  increase_plot_area=1.1)
p1_ancient_WGCNA
```


### Total

```{r}
p1_total_WGCNA_TOP25 <- plot_enrichResult_table_in_network(
  enrich_res_table_tmp = total_WGCNA,
  title = "total_WGCNA",
  output_path=paste(save_dir_tmp, "/test_representation_cluster_total_WGCNA_with_GO_name.pdf", sep = ""),
  width_output=12,
  height_output=12,
  resolution_clustering=1.2,
  threshold_percent_edges = 0.2,
  force_attraction_cluster=4,
  force_repulsion_name=10,
  seed_tmp=123,
  plot_edge_number=F,
  plot_GO_name = T,
  increase_plot_area=1.1)
p1_total_WGCNA
```

```{r}
p1_total_WGCNA_TOP25 <- plot_enrichResult_table_in_network(
  enrich_res_table_tmp = total_WGCNA %>% arrange(p.adjust) %>% head(25),
  title = "total_WGCNA_TOP25",
  output_path=paste(save_dir_tmp, "/test_representation_cluster_total_WGCNA_TOP25_with_GO_name.pdf", sep = ""),
  width_output=12,
  height_output=12,
  resolution_clustering=1.2,
  threshold_percent_edges = 0.2,
  force_attraction_cluster=4,
  force_repulsion_name=10,
  seed_tmp=123,
  plot_edge_number=F,
  plot_GO_name = T,
  increase_plot_area=1.1)
p1_total_WGCNA_TOP25
```

```{r}
p1_total_DE <- plot_enrichResult_table_in_network(
  enrich_res_table_tmp = total_DE,
  output_path=paste(save_dir_tmp, "/test_representation_cluster_total_DE_with_GO_name.pdf", sep = ""),
  title = "total_DE",
  width_output=6,
  height_output=6,
  resolution_clustering=1,
  threshold_percent_edges = 0.1,
  force_attraction_cluster=4,
  force_repulsion_name=10,
  seed_tmp=123,
  plot_edge_number=F,
  plot_GO_name = T,
  increase_plot_area=1.1)
p1_total_DE
```





### Recent

```{r}
p1_recent_DE <- plot_enrichResult_table_in_network(
  enrich_res_table_tmp = recent_DE,
  output_path=paste(save_dir_tmp, "/test_representation_cluster_recent_DE_with_GO_name.pdf", sep = ""),
  width_output=6,
  height_output=6,
  title = "recent_DE",
  resolution_clustering=1,
   threshold_percent_edges = 0.1,
  force_attraction_cluster=4,
  force_repulsion_name=10,
  seed_tmp=123,
  plot_edge_number=F,
  plot_GO_name = T,
  increase_plot_area=1.1)
p1_recent_DE
```

```{r}
p1_recent_WGCNA <- plot_enrichResult_table_in_network(
  enrich_res_table_tmp = recent_WGCNA,
  output_path=paste(save_dir_tmp, "/test_representation_cluster_recent_WGCNA_with_GO_name.pdf", sep = ""),
  title = "recent_WGCNA",
  width_output=15,
  height_output=15,
  resolution_clustering=1.2,
  threshold_percent_edges = 0.2,
  force_attraction_cluster=10,
  force_repulsion_name=40,
  seed_tmp=123,
  plot_edge_number=F,
  plot_GO_name = T,
  increase_plot_area=1.1)
p1_recent_WGCNA
```



```{r}
p1_recent_WGCNA_TOP25 <- plot_enrichResult_table_in_network(
  enrich_res_table_tmp = recent_WGCNA%>% arrange(p.adjust) %>% head(25),
  output_path=paste(save_dir_tmp, "/test_representation_cluster_recent_WGCNA_TOP25_with_GO_name.pdf", sep = ""),
  title = "recent_WGCNA_TOP25",
  width_output=6,
  height_output=6,
  resolution_clustering=1,
  threshold_percent_edges = 0.2,
  force_attraction_cluster=4,
  force_repulsion_name=40,
  seed_tmp=123,
  plot_edge_number=F,
  plot_GO_name = T,
  increase_plot_area=1.1)
p1_recent_WGCNA_TOP25
```

### Murinae

```{r}
p1_murinae_DE <- plot_enrichResult_table_in_network(
  enrich_res_table_tmp = murinae_DE,
  output_path=paste(save_dir_tmp, "/test_representation_cluster_murinae_DE_with_GO_name.pdf", sep = ""),
  width_output=8,
  height_output=8,
  title = "murinae_DE",
  resolution_clustering=1,
  threshold_percent_edges = 0.1,
  force_attraction_cluster=4,
  force_repulsion_name=10,
  seed_tmp=123,
  plot_edge_number=F,
  plot_GO_name = T,
  increase_plot_area=1.1)
p1_murinae_DE
```

```{r fig.height=6, fig.width=6}
p1_murinae_WGCNA_TOP25 <- plot_enrichResult_table_in_network(
enrich_res_table_tmp = murinae_WGCNA %>% arrange(p.adjust) %>% head(25),
output_path=paste(save_dir_tmp, "/test_representation_cluster_murinae_WGCNA_TOP25_with_GO_name.pdf", sep = ""),
title = "murinae_WGCNA_TOP25",
width_output=6,
height_output=6,
resolution_clustering=1,
threshold_percent_edges = 0.1,
force_attraction_cluster=2,
force_repulsion_name=40,
seed_tmp=123,
plot_edge_number=F,
plot_GO_name = T,
increase_plot_area=1.1)
p1_murinae_WGCNA_TOP25
```


```{r fig.height=6, fig.width=6}
p1_murinae_WGCNA <- plot_enrichResult_table_in_network(
enrich_res_table_tmp = murinae_WGCNA,
output_path=paste(save_dir_tmp, "/test_representation_cluster_murinae_WGCNA_with_GO_name.pdf", sep = ""),
title = "murinae_WGCNA",
width_output=18,
height_output=18,
resolution_clustering=1.5,
threshold_percent_edges = 0.2,
force_attraction_cluster=10,
force_repulsion_name=40,
seed_tmp=123,
plot_edge_number=F,
plot_GO_name = T,
increase_plot_area=1.1)
p1_murinae_WGCNA
```

```{r}
list_p1 <- lapply(c(p1_total_WGCNA,
                  p1_total_DE,
                  p1_ancient_WGCNA,
                  p1_ancient_DE,
                  p1_recent_WGCNA,
                  p1_recent_DE,
                  p1_murinae_WGCNA,
                  p1_murinae_DE),
                  function(p){p[[1]]})
```





# sessionInfo

```{r}
sessionInfo()
```

