11.modularity_analysis_shuffled.R

# includes -----------
library(tidyverse)
library(dplyr)
library(magrittr)
library(igraph)
library(reshape2)
library(cowplot)#?
library(vegan)
library(data.table)
source('functions.R')

# Get observed modularity values ----
# read observed for comparing later:
modules_obs <- read_csv('local_output/farm_modules_pos_30_U.csv')
mod_summary_obs <- read_csv('local_output/farm_modules_pos_30_summary.csv', 
                            col_names = c('net', 'call', 'L', 'top_modules', 'time_stamp'))
# get the latest run
mod_summary_obs <- mod_summary_obs[which.max(mod_summary_obs$time_stamp),] 
num_modules_obs <- mod_summary_obs$top_modules
L_obs <- mod_summary_obs$L

# read shuffled *summary* for comparing later:
parent.folder <- "HPC/shuffled/shuffle_farm_r0_30_500_jac_intra"

# General stats
summ <- read_csv(paste(parent.folder,'/farm_modulation_summary_pf_unif.csv',sep=''), 
                 col_names = c('e_id','JOB','call','L','num_modules'))
summ %<>% slice(tail(row_number(), 1000)) %>% filter(str_detect(call, '-2')) # We want only runs with *no* multi-level analysis
ggplot(summ, aes(L))+
  geom_histogram()+
  geom_vline(xintercept = L_obs, linetype = "dashed")
pvalue <- sum(summ$L<L_obs)/500

ggplot(summ, aes(num_modules))+
  geom_histogram()+
  geom_vline(xintercept = num_modules_obs, linetype = "dashed")

# percent of shuffled networks with 1 big module:
100*sum(summ$num_modules == 1)/nrow(summ)

## Read shuffled network themselves -------------------------------
# Sub-folders
sub.folders <- list.dirs(parent.folder, recursive=TRUE)[-1]

# get model nums (because there was NA in the summ for some reason)
module_nums <- NULL

# read all files 
modules_shuffled <- NULL
for (s in sub.folders) {
  print(s)
  modules_run <- list.files(path = s , pattern = '_farm_modules_pf_unif.csv', recursive = T,full.names = T)
  modules_run <- fread(modules_run)
  modules_run$id <- str_split_fixed(s, pattern = '/', n = 4)[4]
  modules_shuffled <- rbind(modules_shuffled,modules_run)
  module_nums <- rbind(module_nums, data.table(origin=s, modules=max(modules_run$module)))
}
modules_shuffled <- as_tibble(modules_shuffled)

write_csv(modules_shuffled, 'local_output/farm_modules_shuffled_500.csv')


# percent of shuffled networks with 1 big module as read from the modules files:
100* (module_nums$modules == 1)/nrow(module_nums)


a2 <- ggplot(summ, aes(L))+
  geom_histogram()+
  geom_vline(xintercept = L_obs, linetype = "dashed", color= 'red') + 
  paper_figs_theme

b2 <- ggplot(module_nums, aes(modules))+
  geom_histogram()+ xlab("Number of modules") +
  geom_vline(xintercept = num_modules_obs, linetype = "dashed", color= 'red') + 
  paper_figs_theme

pdf('local_output/figures/SI_modularity_obs_shuff.pdf',10,6)
plot_grid(a2, b2, nrow = 1, ncol = 2, labels = c('(A)','(B)'),vjust = 1.1, hjust = 0)
dev.off()


# explore shuffled networks ----
farm_modules_shuffled <- read_csv('local_output/farm_modules_shuffled_500.csv')

# number of modules in each permutation
farm_modules_shuffled %>%
  group_by(id) %>%
  summarise(module_sum=n_distinct(module)) %>%
  ggplot(aes(id,module_sum)) +
  geom_col() +
  theme_bw() +
  labs(x='', y='', title='Number of modules in each permutation')+
  theme(axis.text.x = element_text(angle = 90, size=6, color='black'))

n_distinct(farm_modules_shuffled$module)

# number of layers in modules
farm_modules_shuffled %>%
  group_by(id, module) %>%
  summarise(number_of_farms=n_distinct(layer_name)) %>%
  ggplot(aes(number_of_farms)) +
  geom_histogram(color='white')+
  labs(x='Number of farms', title='Number of farms in modules')+
  theme(axis.text.x = element_text(size=10, color='black'))+
  theme_bw()

# nodes in modules
farm_modules_shuffled %>%
  group_by(module) %>%
  summarise(physical_nodes=n_distinct(node_name)) %>%
  arrange(desc(physical_nodes)) %>%
  ggplot(aes(physical_nodes)) +
  geom_histogram(color='white') +
  theme_bw()+
  labs(x='Number of physical nodes', title='Nodes in modules (shuffled networks)')+
  theme(axis.text = element_text(size=15, color='black'), axis.title = element_text(size=15, color='black'), title = element_text(size = 16))

# number of physical nodes in modules
farm_modules_shuffled %>%
  group_by(id, module) %>%
  summarise(physical_nodes=n_distinct(node_name)) %>%
  arrange(desc(physical_nodes)) %>%
  ggplot(aes(physical_nodes)) +
  geom_histogram(color='white')+ scale_y_log10() +
  theme_bw()+
  labs(x='Number of physical nodes', title='Number of physical nodes in modules (shuffled between)')+
  theme(axis.text = element_text(size=12, color='black'), axis.title = element_text(size=15, color='black'), title = element_text(size = 16))

# looking at the biggest module----
# observed
farm_modules_pos <- read_csv('local_output/farm_modules_pos_30_U.csv')
observed_biggest_module <- farm_modules_pos %>%
  group_by(module) %>%
  mutate(layers_in_modules=n_distinct(layer_name)) %>%
  arrange(desc(layers_in_modules)) %>%
  filter(layers_in_modules==max(layers_in_modules)) %>%
  mutate(nodes_in_layers=n_distinct(node_name)) %>%
  arrange(desc(nodes_in_layers)) %>%
  filter(module == 1) %>%
  mutate(id=0) %>%
  select(c(id,module, node_name, layer_name=short_name, layers_in_modules, nodes_in_layers)) %>%
  transform(id = as.character(id)) %>%
  tibble()
# shuffled
shuffled_biggest_module <- farm_modules_shuffled %>%
  select(-c(node_id,layer_id, flow)) %>%
  group_by(id,module) %>%
  mutate(layers_in_modules=n_distinct(layer_name)) %>%
  arrange(desc(layers_in_modules)) %>%
  filter(layers_in_modules==7) %>%
  group_by(id,module,layer_name) %>%
  mutate(nodes_in_layers=n_distinct(node_name)) %>%
  arrange(desc(nodes_in_layers)) %>%
  filter(module == 1) %>%
  transform(module = as.integer(module)) %>%
  tibble()

# combine observed and shuffled
observed_biggest_module$group <- 'obs'
shuffled_biggest_module$group <- 'shuff'
biggest_module_combined <- rbind(observed_biggest_module,shuffled_biggest_module)

obs <- biggest_module_combined %>% 
  group_by(group,id,layer_name) %>% 
  summarise(n=n_distinct(node_name)) %>% 
  filter(group=='obs')

png(filename = 'local_output/figures/observed_shuffled_biggest_module.png', width = 1600, height = 900, res = 300)
biggest_module_combined %>% 
  group_by(group,id,layer_name) %>% 
  summarise(n=n_distinct(node_name)) %>% 
  filter(group=='shuff') %>% 
  ggplot(aes(n))+geom_histogram()+
  facet_wrap(~layer_name)+ 
  theme_bw()+
  geom_vline(data=obs, aes(xintercept = n), color= 'red')+
  theme(axis.text = element_text(size=8 ,color='black'), axis.title = element_text(size=18 ,color='black'),title = element_text(size=16 ,color='black'))+
  labs(x = 'number of nodes',y='count', title = 'Observed vs Shuffled (large module)')
dev.off()

# intercept
biggest_module_combined %>%
  group_by(layer_name,group,id) %>%
  summarise(number_of_nodes=n_distinct(node_name)) %>%
  filter(layer_name=='UK1')

# combine all modules
farm_modules_observed <- farm_modules_pos %>%
  mutate(id='000') %>%
  select(c(id,module, node_name, layer_name=short_name)) %>%
  transform(id = as.character(id)) %>%
  tibble()

farm_modules_shuffled <- farm_modules_shuffled %>%
  select(-c(node_id,layer_id)) %>%
  transform(module = as.integer(module)) %>%
  tibble()

# combine observed and shuffled
farm_modules_observed$group <- 'obs'
farm_modules_shuffled$group <- 'shuff'
farm_modules_combined <- rbind(farm_modules_observed,
                               farm_modules_shuffled %>% select(id, module, node_name, layer_name, group))
  
# comparison between observed and shuffled - percentage in each farm
farm_modules_filtered <- farm_modules_combined %>%
  group_by(group,id,node_name) %>%
  mutate(number_of_farms=n_distinct(layer_name)) %>%
  filter(number_of_farms==7)

biggest_module_filtered <- farm_modules_filtered %>%
  filter(module==1)

farm_modules_filtered %<>%
  group_by(group,id,layer_name) %>%
  summarise(number_of_nodes=n_distinct(node_name))

biggest_module_filtered %<>%
  group_by(group,id,layer_name) %>%
  summarise(number_of_nodes_big_module=n_distinct(node_name))

ASVs_big_module_farm <- left_join(farm_modules_filtered,biggest_module_filtered) 

ASVs_big_module_farm %<>%
  mutate(percent = 100 * number_of_nodes_big_module/number_of_nodes ) %>%
  select(-c(number_of_nodes,number_of_nodes_big_module)) %>%
  group_by(layer_name,group) %>%
  summarise(percent_mean=mean(percent))

png(filename = 'local_output/figures/nodes_occur_in_all_farms.png', width = 1600, height = 900, res = 300)
ASVs_big_module_farm %>%
  ggplot(aes(x=layer_name, y=percent_mean, fill=group)) +
  geom_bar(stat="identity", position="dodge")+
  theme_bw()+
  labs(x = 'Layer',y='percent of nodes in the largest module', title = 'Percent of nodes in the largest module')+
  theme(axis.text = element_text(size=10 ,color='black'), axis.title = element_text(size=12 ,color='black'), title = element_text(size=10 ,color='black'), ) 
dev.off()

# percent of nodes in the big module
farm_modules_pos %>%
  mutate(number_of_nodes=n_distinct(node_id)) %>%
  filter(module==1) %>%
  mutate(number_of_nodes_big_m=n_distinct(node_id)) %>%
  mutate(percent = 100 * number_of_nodes_big_m/number_of_nodes) %>%
  distinct(percent)

# filter the smallest modules----
# for the observed network
farm_modules_pos_03 <- read_csv('local_output/farm_modules_pos_30_U.csv')
farm_modules_pos_03 %>%
  group_by(module) %>%
  summarise(n=n_distinct(node_id)) %>%
  filter(n>1) 

# for the shuffled networks
farm_modules_shuff_filtered <- farm_modules_shuffled %>%
  group_by(id,module) %>% 
  summarise(n=n_distinct(node_id)) %>%
  filter(n>1) 

png(filename = 'local_output/figures/observed_shuffled_hist_filter.png', width = 1600, height = 900, res = 300)
farm_modules_shuff_filtered %>%
  group_by(id) %>%
  summarise(module_sum=n_distinct(module)) %>%
  ggplot(aes(module_sum)) +
  geom_histogram() +
  theme_bw() +
  geom_vline(aes(xintercept=8), color= 'red') + 
  geom_text(aes(x=7.6,label='Observed', y=450, angle=90)) +
  theme(axis.text.x = element_text(size=10, color='black'),
        axis.text.y = element_text(size = 10, color='black'))+
  labs(y='count', title = 'Observed vs Shuffled (filter)')
dev.off()

# pairwise similarity big module-----
# observed:
presence_absence_big_module <- observed_biggest_module %>%
  distinct(layer_name, node_name) %>%
  mutate(present=1) %>% 
  spread(node_name, present, fill = 0)  %>%
  column_to_rownames('layer_name')
presence_absence_big_module <- as.matrix(presence_absence_big_module)
dim(presence_absence_big_module)
beta_div_big_module <- 1-as.matrix(vegdist(presence_absence_big_module, "jaccard"))
# Heatmap 
beta_div_big_module_m <- melt(beta_div_big_module)

ggplot(beta_div_big_module_m, aes(x = Var1, y = Var2, fill = value, label=round(value,2))) +
  geom_tile() + 
  geom_text()+
  theme_bw()+
  scale_fill_gradient(high = "blue", low = "light blue") +
  labs(x = '',y='', title = 'Shared nodes between farms (big module)')+
  theme(axis.text = element_text(size=14, color='black'), title = element_text(size = 14, color = 'black'))

beta_div_big_module_m %<>%
  tibble() %>%
  filter(value!=1) %>%
  distinct(value, .keep_all = TRUE) %>%
  mutate(id='000') %>%
  unite(farm_pairs,Var1:Var2, sep = '-')

# shuffled farms:
shuffled_biggest_module_sumry <- shuffled_biggest_module %>%
  distinct(id,layer_name, node_name) %>%
  mutate(present=1) %>%
  arrange(id)

# first we create a tibble includes all the farm pairs:
# a loop for all ids
beta_div_big_module_shuffled_final <- NULL
for (i in unique(shuffled_biggest_module_sumry$id)) { 
  print(i)
  presence_absence_big_module_shuffled <- shuffled_biggest_module_sumry %>%
    filter(id==i) %>%
    spread(node_name, present, fill = 0)  %>%
    select(-id) %>%
    column_to_rownames('layer_name')
  presence_absence_big_module_shuffled <- as.matrix(presence_absence_big_module_shuffled)
  beta_div_big_module_shuffled <- 1-as.matrix(vegdist(presence_absence_big_module_shuffled, "jaccard"))
  beta_div_big_module_shuffled_m <- melt(beta_div_big_module_shuffled)
  beta_div_big_module_shuffled_m %<>%
    tibble() %>%
    filter(value!=1) %>%
    distinct(value, .keep_all = TRUE) %>%
    mutate(id=i) %>%
    unite(farm_pairs,Var1:Var2, sep = '-')
  beta_div_big_module_shuffled_final <- rbind(beta_div_big_module_shuffled_final,beta_div_big_module_shuffled_m)
}

beta_div_big_module_m$group <- 'obs'
beta_div_big_module_shuffled_final$group <- 'shuff'
beta_div_big_module_combined <- rbind(beta_div_big_module_m,beta_div_big_module_shuffled_final)

beta_div_big_module_combined %>%
  ggplot(aes(value, fill=group)) +
  geom_histogram(binwidth = 2, position = 'dodge') +
  theme_bw() +
  facet_wrap(~farm_pairs)+
  labs(x = 'nodes similarity',y='count', title = 'Observed vs Shuffled nodes similarity')

# for each pair of farms
beta_div_big_module_shuffled_final %>%
  filter(farm_pairs=='SE1-FI1') %>%
  ggplot(aes(value)) +
  geom_histogram(color='white',binwidth = 0.008, position = 'dodge') +
  theme_bw() +
  geom_vline(aes(xintercept=0.0326), color= 'red')+
  labs(x = 'nodes similarity',y='count', title = 'Observed vs Shuffled nodes similarity SE1-FI1')


# analysis with shuffled between farms-----
# combine for biggest module with the between shuffled
observed_biggest_module$group <- 'obs'
shuffled_biggest_module$group <- 'shuff'
biggest_module_combined_all <- rbind(observed_biggest_module,shuffled_biggest_module)

# how many ASVs in each layer
biggest_module_combined_all %<>%
  group_by(layer_name,group) %>%
  summarise(number_of_ASVs_biggest_module=n_distinct(node_name))

# combine all
farm_modules_observed <- farm_modules_pos %>%
  mutate(id='000') %>%
  select(c(id,module, node_name, layer_name=short_name)) %>%
  transform(id = as.character(id)) %>%
  tibble()

farm_modules_shuffled_new <- farm_modules_shuffled %>%
  select(-c(node_id, layer_id, flow)) %>%
  transform(module = as.integer(module)) %>%
  tibble()

# combine observed and shuffled
farm_modules_observed$group <- 'obs'
farm_modules_shuffled_new$group <- 'shuff'
farm_modules_combined <- rbind(farm_modules_observed,farm_modules_shuffled_new)

# how many ASVs in each layer
farm_modules_combined_ASVs <- farm_modules_combined %>%
  group_by(layer_name,group) %>%
  summarise(number_of_ASVs=n_distinct(node_name))  

# join farm_modules_combined and biggest_module_combined_all
farm_modules_comb <- left_join(biggest_module_combined_all,farm_modules_combined_ASVs)
farm_modules_comb %<>%
  mutate(relative_ASVs=number_of_ASVs_biggest_module/number_of_ASVs)

farm_modules_comb %>%
  ggplot(aes(x=layer_name,y=relative_ASVs, fill=group)) +
  geom_bar(stat="identity", position="dodge")+
  theme_bw()+
  theme(axis.text = element_text(size=10, color='black')) 

# how many ASVs that occurr in all 7 farms, are included in the big module 
# comparison between observed and shuffled
farm_modules_combined_filtered <- farm_modules_combined %>%
  group_by(group,id,node_name) %>%
  mutate(number_of_farms=n_distinct(layer_name)) %>%
  filter(number_of_farms==7)

biggest_module_combined_filterd <- farm_modules_combined_filtered %>%
  filter(module==1)

farm_modules_combined_filtered %<>%
  group_by(group,id) %>%
  summarise(number_of_nodes=n_distinct(node_name))

biggest_module_combined_filterd %<>%
  group_by(group,id) %>%
  summarise(number_of_nodes_big_module=n_distinct(node_name))

ASVs_big_module <- left_join(farm_modules_combined_filtered,biggest_module_combined_filterd) 

ASVs_big_module %<>% melt(id.vars= c("group", "id"))   # converting wide data to long

png(filename = 'output/figures/ASVs_occur_in_all_farms.png', width = 1600, height = 900, res = 300)
ASVs_big_module %>%
  ggplot(aes(x=group, y=value, fill=variable)) +
  geom_bar(stat="identity", position="dodge")+
  theme_bw() +
  labs(x = '',y='number of nodes', title = 'ASVs that occur in all 7 farms, largest module vs all modules')+
  theme(axis.text = element_text(size=9 ,color='black'), axis.title = element_text(size=12 ,color='black'), title = element_text(size=10 ,color='black'), )
dev.off()