.Rhistory

nodes_in_tree <- tree$tip.label
# print(setdiff(nodes_in_modules, nodes_in_tree)) # In modules but not in tree
# print(setdiff(nodes_in_tree, nodes_in_modules)) # In tree but not in modules
# Overlapping nodes:
overlapping <- intersect(nodes_in_tree, nodes_in_modules)
# Observed modules
M_obs <- module_object$modules %>%
filter(str_starts(node_name, node_start_letter)) %>%
mutate(node_name=str_replace_all(node_name, pattern = '\\.', '')) %>%
filter(node_name %in% overlapping) %>%
rename(m=module_level1) %>%
select(node_name, m)
#Mean PDistance between hosts within modules
D_obs <- M_obs %>%
inner_join(D, by=c('node_name'='from')) %>% # join PD distances
rename(d=weight) %>%
arrange(m, node_name) %>%
group_by(m) %>% # Per module
filter(to %in% node_name) %>% #Host pairs within a module
summarise(d_mean=mean(d), mod_size=n())
D_obs_mean <- mean(D_obs$d_mean)
# print('Observed network:')
print(D_obs)
#Shuffle to create permuted modules of the same size,
#and recalculate the meand PD within modules. The shuffling permutes the ID of the strains.
D_perm <- NULL
nperm <- 500
for (i in 1:nperm){
# print(i)
D_perm %<>% bind_rows(
M_obs %>%
mutate(node_name=sample(node_name, replace = F)) %>%
inner_join(D, by=c('node_name'='from')) %>% # join PD distances
rename(d=weight) %>%
arrange(m, node_name) %>%
group_by(m) %>% # Per module
filter(to %in% node_name) %>% #Host pairs within a module
summarise(d_mean=mean(d)) %>% # Calculate mean PD within modules
mutate(run=i)
)
}
# Null hypothesis is that the permuted distance is smaller than the observed for
# each module (i.e., no signal). If we reject this hypothesis then there is
# phylogenetic signal because the observed PD beteween hosts within each module
# would be smaller than expected by chance (closely related hosts share a module).
# Plot the means
plt_across_modules <-
D_perm %>% group_by(run) %>%
summarise(D_perm_mean=mean(d_mean)) %>%
ggplot(aes(x=D_perm_mean))+geom_histogram()+geom_vline(xintercept = D_obs_mean)
result_across_moduels <-
D_perm %>% group_by(run) %>%
summarise(D_perm=mean(d_mean)) %>%
mutate(test=D_perm<D_obs_mean) %>%
summarise(pvalue=sum(test)/nperm) %>%
mutate(res=ifelse(pvalue<0.05,'Signal','No signal'))
# This can also be tested per module
plt_within_modules <-
D_perm %>%
full_join(D_obs, by='m') %>%
rename(d_perm=d_mean.x, d_obs=d_mean.y) %>%
ggplot(aes(x=d_perm))+
geom_histogram()+
facet_wrap(~m)+
geom_vline(data = D_obs, aes(xintercept = d_mean))
result_within_moduels <-
D_perm %>%
full_join(D_obs, by='m') %>%
rename(d_perm=d_mean.x, d_obs=d_mean.y) %>%
mutate(test=d_perm<d_obs) %>%
group_by(m) %>%
summarise(pvalue=sum(test)/nperm) %>%
mutate(Signif=ifelse(pvalue<0.05,'Signal','No signal'),
Signif_Bonferroni=ifelse(pvalue<0.05/nrow(D_obs),'Signal','No signal')) # Need to divide by number of modules for Bonferroni correction
out <- list(D_obs=D_obs,
D_obs_mean=D_obs_mean,
plt_across_modules=plt_across_modules,
plt_within_modules=plt_within_modules,
result_across_moduels=result_across_moduels,
result_within_moduels=result_within_moduels,
nodes_in_modules=nodes_in_modules,
nodes_in_tree=nodes_in_tree,
overlapping=overlapping)
return(out)
}
i=1
hr <- VDRs$end[i]
edges <- all_networks[[which(hr_seq==hr)]]$infection_list
x <- create_monolayer_object(edges, directed = F, bipartite = T)
test <- run_infomap_monolayer_nonrandom(x, infomap_executable = 'Infomap', flow_model = 'undirected', silent = T, trials = 100, two_level = T, seed = 123, signif = T, shuff_method = 'r00', nsim = 100)
infection_modularity %<>% bind_rows(tibble(hr=hr,
pvalue=test$pvalue,
n_hosts=ncol(x$mat),
n_spacer=nrow(x$mat),
n_interactions=sum(x$mat>0),
n_modules=max(test$modules$module_level1)))
test <- run_infomap_monolayer(x, infomap_executable = 'Infomap', flow_model = 'undirected', silent = T, trials = 100, two_level = T, seed = 123, signif = T, shuff_method = 'r00', nsim = 100)
infection_modularity %<>% bind_rows(tibble(hr=hr,
pvalue=test$pvalue,
n_hosts=ncol(x$mat),
n_spacer=nrow(x$mat),
n_interactions=sum(x$mat>0),
n_modules=max(test$modules$module_level1)))
plot_modular_matrix(test)
tibble(L_sim=test$L_sim) %>%
ggplot(aes(L_sim))+
geom_histogram(fill='plum')+
geom_vline(xintercept = test$L, linetype='dashed')+
theme_bw()+
labs(x='Map equation L', y='Count')+
theme(legend.position='none', axis.text = element_text(size=20), axis.title = element_text(size=20))
tree_data <- read_delim(paste(base_name,'_Phage-TREE.txt',sep=''), delim = '\t',col_names = c("Recordtime","id","parent_id","creation_time"))
tree_data$id <- paste('V_',str_pad(tree_data$id, 4, 'left', '0'),sep='')
tree_data$parent_id <- paste('V_',str_pad(tree_data$parent_id, 4, 'left', '0'),sep='')
tree_data %<>% left_join(virus_dynamics_list %>% select(V_ID,death), by=c('id'='V_ID'))
tree_data[1,3] <- NA
tree_data$death[is.na(tree_data$death)]<-tree_data$creation_time[is.na(tree_data$death)]+1
# These data frames contain the extinction/mutation events and persistence of viruses/bacteria
virus_dynamics_list <- virus_data %>%
filter(timesOfRecord<=max(hr_seq)) %>%
select(timesOfRecord, label) %>%
rename(V_ID=label) %>%
mutate(V_ID=paste('V_',str_pad(V_ID, 4, 'left', '0'),sep=''), w=1) %>%
arrange(timesOfRecord, V_ID)
virus_dynamics_list %<>% group_by(V_ID) %>% summarise(birth=first(timesOfRecord), death=last(timesOfRecord)) %>%
mutate(persistence=death-birth+1)
virus_dynamics_list %<>% filter(birth>=min(hr_seq), death<=max(hr_seq))
bacteria_dynamics_list <-
bacteria_data %>%
filter(timesOfRecord<=max(hr_seq)) %>%
select(timesOfRecord, label) %>%
rename(B_ID=label) %>%
mutate(B_ID=paste('B_',str_pad(B_ID, 4, 'left', '0'),sep=''), w=1) %>%
arrange(timesOfRecord, B_ID)
bacteria_dynamics_list %<>% group_by(B_ID) %>% summarise(birth=first(timesOfRecord), death=last(timesOfRecord)) %>%
mutate(persistence=death-birth+1)
bacteria_dynamics_list %<>% filter(birth>=min(hr_seq), death<=max(hr_seq))
tree_data <- read_delim(paste(base_name,'_Phage-TREE.txt',sep=''), delim = '\t',col_names = c("Recordtime","id","parent_id","creation_time"))
tree_data$id <- paste('V_',str_pad(tree_data$id, 4, 'left', '0'),sep='')
tree_data$parent_id <- paste('V_',str_pad(tree_data$parent_id, 4, 'left', '0'),sep='')
tree_data %<>% left_join(virus_dynamics_list %>% select(V_ID,death), by=c('id'='V_ID'))
tree_data[1,3] <- NA
tree_data$death[is.na(tree_data$death)]<-tree_data$creation_time[is.na(tree_data$death)]+1
tree <- nodes_dataframe_to_newick(tree_data)
writeLines(tree, 'tree.nwk')
tree_viruses <- treeio::read.tree('tree.nwk')
plt_viruses_tree <-
standard_plot(
ggtree::ggtree(tree_viruses, ladderize = T) +
ggtree::theme_tree2()
)
tree_data <- read_delim(paste(base_name,'_Bacteria-TREE.txt',sep=''), delim = '\t',col_names = c("Recordtime","id","parent_id","creation_time"))
tree_data$id <- paste('B_',str_pad(tree_data$id, 4, 'left', '0'),sep='')
tree_data$parent_id <- paste('B_',str_pad(tree_data$parent_id, 4, 'left', '0'),sep='')
tree_data %<>% left_join(bacteria_dynamics_list %>% select(B_ID,death), by=c('id'='B_ID'))
tree_data[1,3] <- NA
tree_data$death[is.na(tree_data$death)]<-tree_data$creation_time[is.na(tree_data$death)]+1
tree <- nodes_dataframe_to_newick(tree_data)
writeLines(tree, 'tree_bacteria.nwk')
tree_bacteria <- treeio::read.tree('tree_bacteria.nwk')
x <- test_PD_modules(bac)
x <- test_PD_modules(tree_bacteria, module_object = test, node_start_letter = 'B')
tree = tree_bacteria
module_object = test
node_start_letter = 'B'
# Phylogenetic signal analysis
D <- ape::cophenetic.phylo(tree) # Phyloegentic distance
D <- matrix_to_list_unipartite(D, directed = T) # Use directed to make sure that the from column has all the nodes (need it for joining later)
D <- D$edge_list
# Difference between tree and matrix
nodes_in_modules <- module_object$modules %>%
filter(str_starts(node_name, node_start_letter)) %>%
distinct(node_name) %>%
mutate(node_name=str_replace_all(node_name, pattern = '\\.', ''))
nodes_in_modules <- nodes_in_modules$node_name
nodes_in_tree <- tree$tip.label
# print(setdiff(nodes_in_modules, nodes_in_tree)) # In modules but not in tree
# print(setdiff(nodes_in_tree, nodes_in_modules)) # In tree but not in modules
# Overlapping nodes:
overlapping <- intersect(nodes_in_tree, nodes_in_modules)
print(setdiff(nodes_in_modules, nodes_in_tree)) # In modules but not in tree
print(setdiff(nodes_in_tree, nodes_in_modules)) # In tree but not in modules
overlapping
# Observed modules
M_obs <- module_object$modules %>%
filter(str_starts(node_name, node_start_letter)) %>%
mutate(node_name=str_replace_all(node_name, pattern = '\\.', '')) %>%
filter(node_name %in% overlapping) %>%
rename(m=module_level1) %>%
select(node_name, m)
#Mean PDistance between hosts within modules
D_obs <- M_obs %>%
inner_join(D, by=c('node_name'='from')) %>% # join PD distances
rename(d=weight) %>%
arrange(m, node_name) %>%
group_by(m) %>% # Per module
filter(to %in% node_name) %>% #Host pairs within a module
summarise(d_mean=mean(d), mod_size=n())
D_obs_mean <- mean(D_obs$d_mean)
# print('Observed network:')
print(D_obs)
#Shuffle to create permuted modules of the same size,
#and recalculate the meand PD within modules. The shuffling permutes the ID of the strains.
D_perm <- NULL
D_obs
i=2
hr <- VDRs$end[i]
edges <- all_networks[[which(hr_seq==hr)]]$infection_list
x <- create_monolayer_object(edges, directed = F, bipartite = T)
test <- run_infomap_monolayer(x, infomap_executable = 'Infomap', flow_model = 'undirected', silent = T, trials = 100, two_level = T, seed = 123, signif = T, shuff_method = 'r00', nsim = 100)
infection_modularity %<>% bind_rows(tibble(hr=hr,
pvalue=test$pvalue,
n_hosts=ncol(x$mat),
n_spacer=nrow(x$mat),
n_interactions=sum(x$mat>0),
n_modules=max(test$modules$module_level1)))
plot_modular_matrix(test)
tibble(L_sim=test$L_sim) %>%
ggplot(aes(L_sim))+
geom_histogram(fill='plum')+
geom_vline(xintercept = test$L, linetype='dashed')+
theme_bw()+
labs(x='Map equation L', y='Count')+
theme(legend.position='none', axis.text = element_text(size=20), axis.title = element_text(size=20))
x <- test_PD_modules(tree = tree_bacteria, module_object = test, node_start_letter = 'B')
pd_results <- test_PD_modules(tree = tree_bacteria, module_object = test, node_start_letter = 'B')
pd_results
pd_results$D_obs
pd_results$D_obs_mean
pd_results$result_across_moduels
pd_results$result_within_moduels
i=3
hr <- VDRs$end[i]
edges <- all_networks[[which(hr_seq==hr)]]$infection_list
x <- create_monolayer_object(edges, directed = F, bipartite = T)
test <- run_infomap_monolayer(x, infomap_executable = 'Infomap', flow_model = 'undirected', silent = T, trials = 100, two_level = T, seed = 123, signif = T, shuff_method = 'r00', nsim = 100)
infection_modularity %<>% bind_rows(tibble(hr=hr,
pvalue=test$pvalue,
n_hosts=ncol(x$mat),
n_spacer=nrow(x$mat),
n_interactions=sum(x$mat>0),
n_modules=max(test$modules$module_level1)))
pd_results <- test_PD_modules(tree = tree_bacteria, module_object = test, node_start_letter = 'B')
pd_results$D_obs
pd_results$result_across_moduels
pd_results$result_within_moduels
# !diagnostics off
# parameters of the model; used by the functions, so need to initialize here
beta=50;phi=10^-7;p=10^-5;q=10^-5;m=0.1;spacer_len=10
# A function to identify if analysis is run on the UChicago's Midway HPC
on_Midway <- function(){ifelse(Sys.getenv('USER')=='pilosofs',T,F)}
# If run with an sbatch pipeline take arguments from there. Otherwise those specified here.
if (length(commandArgs(trailingOnly=TRUE))==0) {
# args <- c('mu5e-7_initialDiffDp1_S50P15_R-13997','5*10^-7','15',F, F)
args <- c('mu1e-7_initialDiffDp1_S10P15_R-12499','1*10^-7','15',F, T)
} else {
args <- commandArgs(trailingOnly=TRUE)
}
base_name <- args[1]
mu <- eval(parse(text = args[2]))
protospacer_len <- eval(parse(text = args[3]))
make_plots <- as.logical(args[4])
complete_analysis <- as.logical(args[5])
if(on_Midway()){system('module load gcc/6.1')}
if(!on_Midway()){setwd(paste('data/',base_name,sep=''))}
setwd("~/GitHub/ecomplab/CRISPR_networks")
# !diagnostics off
# parameters of the model; used by the functions, so need to initialize here
beta=50;phi=10^-7;p=10^-5;q=10^-5;m=0.1;spacer_len=10
# A function to identify if analysis is run on the UChicago's Midway HPC
on_Midway <- function(){ifelse(Sys.getenv('USER')=='pilosofs',T,F)}
# If run with an sbatch pipeline take arguments from there. Otherwise those specified here.
if (length(commandArgs(trailingOnly=TRUE))==0) {
# args <- c('mu5e-7_initialDiffDp1_S50P15_R-13997','5*10^-7','15',F, F)
args <- c('mu1e-7_initialDiffDp1_S10P15_R-12499','1*10^-7','15',F, T)
} else {
args <- commandArgs(trailingOnly=TRUE)
}
base_name <- args[1]
mu <- eval(parse(text = args[2]))
protospacer_len <- eval(parse(text = args[3]))
make_plots <- as.logical(args[4])
complete_analysis <- as.logical(args[5])
if(on_Midway()){system('module load gcc/6.1')}
if(!on_Midway()){setwd(paste('data/',base_name,sep=''))}
if(!on_Midway()){setwd(paste(getwd(),'data/',base_name,sep=''))}
base_name
dom_strains_num
plot_signif(test)
i=1
metaweb <- NULL
for (hr in VDRs$start[i]:VDRs$end[i]){
edges <- all_networks[[which(hr_seq==hr)]]$bacteria_sp_list
if(is.null(edges)){next}
metaweb %<>% bind_rows(edges) %>% distinct(B_ID, SP)
}
metaweb %<>% mutate(w=1)
x <- create_monolayer_object(metaweb, directed = F, bipartite = T)
host_sp_modularity <- run_infomap_monolayer(x, infomap_executable = 'Infomap', flow_model = 'undirected', silent = T,
trials = 100, two_level = T, seed = 123,
signif = F)
PD_results <- test_PD_modules(tree_bacteria, host_sp_modularity, 'B')
# Function to test for phylogenetic signal in modules
test_PD_modules<- function(tree, module_object, node_start_letter){
# Phylogenetic signal analysis
D <- ape::cophenetic.phylo(tree) # Phyloegentic distance
D <- matrix_to_list_unipartite(D, directed = T) # Use directed to make sure that the from column has all the nodes (need it for joining later)
D <- D$edge_list
# Difference between tree and matrix
nodes_in_modules <- module_object$modules %>%
filter(str_starts(node_name, node_start_letter)) %>%
distinct(node_name) %>%
mutate(node_name=str_replace_all(node_name, pattern = '\\.', ''))
nodes_in_modules <- nodes_in_modules$node_name
nodes_in_tree <- tree$tip.label
# print(setdiff(nodes_in_modules, nodes_in_tree)) # In modules but not in tree
# print(setdiff(nodes_in_tree, nodes_in_modules)) # In tree but not in modules
# Overlapping nodes:
overlapping <- intersect(nodes_in_tree, nodes_in_modules)
# Observed modules
M_obs <- module_object$modules %>%
filter(str_starts(node_name, node_start_letter)) %>%
mutate(node_name=str_replace_all(node_name, pattern = '\\.', '')) %>%
filter(node_name %in% overlapping) %>%
rename(m=module_level1) %>%
select(node_name, m)
#Mean PDistance between hosts within modules
D_obs <- M_obs %>%
inner_join(D, by=c('node_name'='from')) %>% # join PD distances
rename(d=weight) %>%
arrange(m, node_name) %>%
group_by(m) %>% # Per module
filter(to %in% node_name) %>% #Host pairs within a module
summarise(d_mean=mean(d), mod_size=n())
D_obs_mean <- mean(D_obs$d_mean)
# print('Observed network:')
# print(D_obs)
#Shuffle to create permuted modules of the same size,
#and recalculate the meand PD within modules. The shuffling permutes the ID of the strains.
D_perm <- NULL
nperm <- 500
for (i in 1:nperm){
# print(i)
D_perm %<>% bind_rows(
M_obs %>%
mutate(node_name=sample(node_name, replace = F)) %>%
inner_join(D, by=c('node_name'='from')) %>% # join PD distances
rename(d=weight) %>%
arrange(m, node_name) %>%
group_by(m) %>% # Per module
filter(to %in% node_name) %>% #Host pairs within a module
summarise(d_mean=mean(d)) %>% # Calculate mean PD within modules
mutate(run=i)
)
}
# Null hypothesis is that the permuted distance is smaller than the observed for
# each module (i.e., no signal). If we reject this hypothesis then there is
# phylogenetic signal because the observed PD beteween hosts within each module
# would be smaller than expected by chance (closely related hosts share a module).
# Plot the means
plt_across_modules <-
D_perm %>% group_by(run) %>%
summarise(D_perm_mean=mean(d_mean)) %>%
ggplot(aes(x=D_perm_mean))+geom_histogram()+geom_vline(xintercept = D_obs_mean)
result_across_moduels <-
D_perm %>% group_by(run) %>%
summarise(D_perm=mean(d_mean)) %>%
mutate(test=D_perm<D_obs_mean) %>%
summarise(pvalue=sum(test)/nperm) %>%
mutate(res=ifelse(pvalue<0.05,'Signal','No signal'))
# This can also be tested per module
plt_within_modules <-
D_perm %>%
full_join(D_obs, by='m') %>%
rename(d_perm=d_mean.x, d_obs=d_mean.y) %>%
ggplot(aes(x=d_perm))+
geom_histogram()+
facet_wrap(~m)+
geom_vline(data = D_obs, aes(xintercept = d_mean))
result_within_moduels <-
D_perm %>%
full_join(D_obs, by='m') %>%
rename(d_perm=d_mean.x, d_obs=d_mean.y) %>%
mutate(test=d_perm<d_obs) %>%
group_by(m) %>%
summarise(pvalue=sum(test)/nperm) %>%
mutate(Signif=ifelse(pvalue<0.05,'Signal','No signal'),
Signif_Bonferroni=ifelse(pvalue<0.05/nrow(D_obs),'Signal','No signal')) # Need to divide by number of modules for Bonferroni correction
out <- list(D_obs=D_obs,
D_obs_mean=D_obs_mean,
plt_across_modules=plt_across_modules,
plt_within_modules=plt_within_modules,
result_across_moduels=result_across_moduels,
result_within_moduels=result_within_moduels,
nodes_in_modules=nodes_in_modules,
nodes_in_tree=nodes_in_tree,
overlapping=overlapping)
return(out)
}
PD_results <- test_PD_modules(tree_bacteria, host_sp_modularity, 'B')
PD_results
PD_results$D_obs
PD_results$D_obs_mean
PD_results$result_within_moduels
print(
PD_results <- test_PD_modules(tree_bacteria, host_sp_modularity, 'B')
PD_results$D_obs
PD_results$result_within_moduels)
print(PD_results$D_obs
PD_results$result_within_moduels)
print(PD_results$D_obs)
print(PD_results$result_within_moduels)
# Significance of modularity of infection networks --------------------------------------
infection_modularity <- NULL
# Significance of modularity of infection networks --------------------------------------
infection_modularity <- NULL
# Test at the end of each VDR
for (i in 1:nrow(VDRs)){
hr <- VDRs$end[i]
edges <- all_networks[[which(hr_seq==hr)]]$infection_list
x <- create_monolayer_object(edges, directed = F, bipartite = T)
test <- run_infomap_monolayer(x, infomap_executable = 'Infomap', flow_model = 'undirected', silent = T, trials = 100, two_level = T, seed = 123, signif = T, shuff_method = 'r00', nsim = 10)
infection_modularity %<>% bind_rows(tibble(hr=hr,
pvalue=test$pvalue,
n_hosts=ncol(x$mat),
n_spacer=nrow(x$mat),
n_interactions=sum(x$mat>0),
n_modules=max(test$modules$module_level1)))
plot_modular_matrix(infection_modularity)
}
plot_modular_matrix(test)
# Significance of modularity of infection networks --------------------------------------
infection_modularity <- NULL
# Test at the end of each VDR
for (i in 1:nrow(VDRs)){
hr <- VDRs$end[i]
edges <- all_networks[[which(hr_seq==hr)]]$infection_list
x <- create_monolayer_object(edges, directed = F, bipartite = T)
test <- run_infomap_monolayer(x, infomap_executable = 'Infomap', flow_model = 'undirected', silent = T, trials = 100, two_level = T, seed = 123, signif = T, shuff_method = 'r00', nsim = 10)
infection_modularity %<>% bind_rows(tibble(hr=hr,
pvalue=test$pvalue,
n_hosts=ncol(x$mat),
n_spacer=nrow(x$mat),
n_interactions=sum(x$mat>0),
n_modules=max(test$modules$module_level1)))
plot_modular_matrix(test)
}
infection_modularity
result_across_moduels
pd_results
pd_results$D_obs
pd_results$result_within_moduels
pd_results$result_within_moduels %>% left_join(pd_results$D_obs)
out <- pd_results$result_within_moduels %>% left_join(pd_results$D_obs)
out$VDR <- i
out
# Phylogenetic signal in infection networks --------------------------------------
phylogenetic_signal_infection <- NULL
# Test at the end of each VDR
for (i in 1:nrow(VDRs)){
hr <- VDRs$end[i]
edges <- all_networks[[which(hr_seq==hr)]]$infection_list
x <- create_monolayer_object(edges, directed = F, bipartite = T)
infection_modularity <- run_infomap_monolayer(x, infomap_executable = 'Infomap', flow_model = 'undirected', silent = T, trials = 100, two_level = T, seed = 123, signif = F)
pd_results <- test_PD_modules(tree = tree_bacteria, module_object = infection_modularity, node_start_letter = 'B')
pd_results$D_obs
out <- pd_results$result_within_moduels %>% left_join(pd_results$D_obs)
out$VDR <- i
phylogenetic_signal_infection <- rbind(phylogenetic_signal_infection, out)
}
phylogenetic_signal_infection
# Phylogenetic signal in host-spacer modules ------------------------------
# Agregate host-spacer networks within each VDR and test for significance of modularity
phylogenetic_signal_hs <- NULL
for (i in 1:nrow(VDRs)){
metaweb <- NULL
for (hr in VDRs$start[i]:VDRs$end[i]){
edges <- all_networks[[which(hr_seq==hr)]]$bacteria_sp_list
if(is.null(edges)){next}
metaweb %<>% bind_rows(edges) %>% distinct(B_ID, SP)
}
metaweb %<>% mutate(w=1)
x <- create_monolayer_object(metaweb, directed = F, bipartite = T)
host_sp_modularity <- run_infomap_monolayer(x, infomap_executable = 'Infomap', flow_model = 'undirected', silent = T,
trials = 100, two_level = T, seed = 123,
signif = F)
pd_results <- test_PD_modules(tree = tree_bacteria, module_object = host_sp_modularity, node_start_letter = 'B')
out <- pd_results$result_within_moduels %>% left_join(pd_results$D_obs)
out$VDR <- i
phylogenetic_signal_hs <- rbind(phylogenetic_signal_hs, out)
}
# Agregate host-spacer networks within each VDR and test for significance of modularity
host_sp_modularity <- NULL
for (i in 1:nrow(VDRs)){
metaweb <- NULL
for (hr in VDRs$start[i]:VDRs$end[i]){
edges <- all_networks[[which(hr_seq==hr)]]$bacteria_sp_list
if(is.null(edges)){next}
metaweb %<>% bind_rows(edges) %>% distinct(B_ID, SP)
}
metaweb %<>% mutate(w=1)
x <- create_monolayer_object(metaweb, directed = F, bipartite = T)
test <- run_infomap_monolayer(x, infomap_executable = 'Infomap', flow_model = 'undirected', silent = T, trials = 100, two_level = T, seed = 123, signif = T, shuff_method = 'r00', nsim = 10)
host_sp_modularity %<>% bind_rows(tibble(hr=hr,
pvalue=test$pvalue,
n_hosts=ncol(x$mat),
n_spacer=nrow(x$mat),
n_interactions=sum(x$mat>0),
n_modules=max(test$modules$module_level1)))
}
record_data(host_sp_modularity)
host_sp_modularity