--- title: "Script_2_Modularity analyses" author: "Jorge" date: "5/18/2022" output: html_document --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) ``` # -------------------------------------------------------------------- # [Title]: Modularity analysis. # Interaction matrices to compute the modularity and node assignment # to modules. # # ```{r load library} library(bipartite) library(ggbipart) library(vcd) library(tibble) library(igraph) options(max.print = 99999) library(bootstrap) library("tcltk") #install.packages("ggrepel") library(ggrepel) library(dplyr) ``` ```{r data loading} mat<-read.csv("~/Documents/GitHub/Juniperus_mutualism_2022/data/visits_mat.csv",header=T,sep=";",dec=".",na.strings="NA") mat<-column_to_rownames(mat, var = "ind") mat<-as.matrix(mat) ``` # Over or infra representation of some stands (plants) in network modules. Sifgificance and results- summary ```{r} # Count number of plants of each stand in each module # Empirical network. res<- computeModules(mat, method="Beckett")# Modularity for observed matrix modlist<- listModuleInformation(res) # Save the observed module # composition data. n<- length(modlist[[2]]) # Number of modules in this run. mm<- as.data.frame(NULL) # Initialize objects. colClasses = c("character", "character") colnames = c("module", "node_stand") mm <- read.table(text = "", colClasses= colClasses, col.names = colnames) for (j in 1:n) { mm1<- cbind(rep(j, length(unlist(modlist[[2]][[j]][[1]]))), substr(unlist(modlist[[2]][[j]][[1]]), start = 1 , stop = 1 )) mm<- rbind(mm, mm1) } colnames(mm)= c("module", "stand.of.node") print(mm) table(mm[,1],mm[,2]) #-------------- #----------------------------------------------------------------------- # Tests # Loop to get the contingency tables. # # Simulated networks. # z<-table(mm[,1],mm[,2]) # Contingency table plot. library("vcd") # plot just a subset of the table assoc(head(z), shade = TRUE, las=3) chi<- chisq.test(z) ; chi; round(chi$residuals, 3) # Residuals of the chi.sq test library(corrplot) corrplot(chi$residuals, is.cor = FALSE) # Contibution in percentage (%) contrib <- 100*chi$residuals^2/chi$statistic round(contrib, 3) corrplot(contrib, is.cor = FALSE) library("graphics") mosaicplot(z, shade = TRUE, las=2, main = "Modules") # library("gplots") # 1. convert the data as a table z <- as.table(as.matrix(z)) # 2. Graph balloonplot(t(z), main ="Modules", xlab ="Stand", ylab="Module #", label = T, show.margins = T) #----------------------------------------------------------------------- ``` # Network plotting # Modularity computation # function loading ```{r m_boot.function} # Function for Bootstrap loop to estimate M for resampled matrices. ------------ # USE WITH: m_boot(mymat, myresamp). myresamp defaults to 99. m_boot <- function (mymat, myresamp= 99) { TIME <- Sys.time() # Resampling plan for bootstrapping n=99 adjacency matrices. m.boot<- NULL # mat is the input matrix for which M is tested # mlike is the observed mean M value (in general, mat@likelihood) # The loop. PUT THIS IN A BOOT LOOP for (i in 1:myresamp) { mat1 <- sample_frac(mymat, 0.8, replace= T) mnmetrics <- computeModules(mat1, method="Beckett", deep= FALSE, deleteOriginalFiles= FALSE, steps= 1000, tolerance= 1e-10, experimental= FALSE, forceLPA= FALSE) m.boot<- rbind(m.boot, mnmetrics@likelihood) } colnames(m.boot)<- c("MBoot") # m.boot saves the M values for bootstrapped # matrices. return(m.boot) # Sys.time() - TIME # } ``` ```{r COL Modularity CI} M_ALL <- computeModules(mat, method="Beckett", deep= FALSE, deleteOriginalFiles= FALSE, steps= 1000, tolerance= 1e-10, experimental= FALSE, forceLPA= FALSE) mat<-as.data.frame(mat) M_ALL@likelihood MObs<- M_ALL@likelihood ALL_Mboot <- as.data.frame(cbind(MObs, # colMeans(m.boot))) colMeans(m.boot<-m_boot(mat, 99)))) rnd.se <- apply(m.boot,2,sd)/sqrt(length(m.boot)) # SE randomized parameters ALL_Mboot<- cbind(ALL_Mboot, rnd.se) ci1 <- MObs-(1.962*ALL_Mboot$rnd.se) ci2 <- MObs+(1.962*ALL_Mboot$rnd.se) ALL_Mboot<- cbind(ALL_Mboot, ci1, ci2) colnames(ALL_Mboot)<- c("MObs", "meanM.boot", "SE", "CI_low", "CI_high") # # Output parameters # meandf # Bootstrap mean of parameter values # SE # Bootstrap SE for parameter values # CI_low # Bootstrap CI (lower) # CI_high # Bootstrap CI (higher) ALL_Mboot ``` ```{r Modularity significance} modintALL <- computeModules(mat, method="Beckett", deep= FALSE, deleteOriginalFiles= FALSE, steps= 1000, tolerance= 1e-10, experimental= FALSE, forceLPA= FALSE) require(bipartite) Msig <- function (mat, mlike) { # mat is the input matrix for which M is tested # mlike is the observed mean M value nulls <- nullmodel(mat, N=100, method=3) modules.nulls <- sapply(nulls, computeModules) like.nulls <- sapply(modules.nulls, function(x) x@likelihood) z <- (mlike - mean(like.nulls))/sd(like.nulls) p <- 2*pnorm(-abs(z)) cat("\n\n","P value for modularity M= ", MObs, "\n", "\n\n", "zeta= ", z, "P= ",format(p, scientific = T),"\n\n") } Msig(mat, MObs) ``` ```{r network plotting} igraph.nev <- graph_from_incidence_matrix(mat, weighted = TRUE) class(igraph.nev) vertex_attr(igraph.nev) edge_attr(igraph.nev) igraph.nev E(igraph.nev) # Species are oredered alphabetically and by trophic level V(igraph.nev) x<-vertex_attr(igraph.nev)$type trophic<-as.factor(ifelse(x=="FALSE","plant","frugivore")) colrs <- c( "yellow", "green") colrs<-adjustcolor(colrs, alpha.f=.8) V(igraph.nev)$color <-colrs[trophic] vertex_attr(igraph.nev) igraph.nev<- delete_vertex_attr(igraph.nev, "color") vertex_attr(igraph.nev) V(igraph.nev)$color <-colrs[trophic] edge_attr(igraph.nev)$weight plot(igraph.nev, edge.width=3*(edge_attr(igraph.nev)$weight)/200, edge.color="black", mark.border=NA) # I use tkplot to create figure 1, in which each node (plant) will occupy its actual location on the landsape. #tkplot(igraph.nev) tplot<- tkplot(igraph.nev) #tkid is the id of the tkplot that will open l <- tkplot.getcoords(tplot) # grab the coordinates from tkplot tk_close(tplot, window.close = F) plot(igraph.nev, layout=l) V(igraph.nev)$size<- 6 plot(igraph.nev,layout=l, edge.width=(edge_attr(igraph.nev)$weight)/400, edge.color="black", mark.border=NA, edge.curved=0.5 ) dev.off() forceNetwork(Links = netnev.w$links, Nodes = netnev.w$nodes, Source = 'source',Target = 'target', NodeID = 'name', Group = 'group', zoom = TRUE, linkDistance = 20,opacity=4,Value ='value',linkColour="grey",charge=-30,legend=T,colourScale =JS("d3.scaleOrdinal(d3.schemeCategory10);")) ``` Estimating standarized -c -z values for visits network ```{r} #Estimating -c -z values cz_animals<-czvalues(res, weighted=TRUE, level="higher") # For animals cz_plants<-czvalues(res, weighted=TRUE, level="lower") # For plants cz_animals<-as.data.frame(cz_animals) cz_plants<-as.data.frame(cz_plants) cz_values<-as.data.frame(rbind(cz_animals,cz_plants)) # Merge plants and animal in a dataset #Estimating percentiles c.crit <- quantile(cz_values$c,probs=c(0.90)) # 0.5969278 z.crit <- quantile(cz_values$z,probs=c(0.90)) # 1.142267 # Traigo datasets para plots #Me traigo dataset-propagule results dispersed_seeds_planttraits<-read.csv("~/Documents/GitHub/Juniperus_mutualism_2022/data/propagule_results.csv",header=T,sep=";",dec=".",na.strings="NA") dispersed_seeds_planttraits<-column_to_rownames(dispersed_seeds_planttraits, var = "ind") # Genero vector para module module<-(dispersed_seeds_planttraits$module) module<-c(c("A","B","B","B","C","C","C","C","C","C","A","C"),module) #Incluyo modulos animales que en el dataset /cz_values/ van delante de las planas modcolors <- c("A"="#b2182b", "B"="#f4a582", "C"="#bababa") # Genero vector para frugivore / plant type_a<-c(rep("frugivore",12)) type_b<-c(rep("plant",105)) type<-c(type_a,type_b) type<-as.data.frame(type) type_color<-c("plant"="#4daf4a" ,"frugivore"="#e41a1c") # Genero vector para area (EN FRUGIVORE LE DEJO FRUGIVORE) col<-c(rep("COL",35)) mar<-c(rep("MAR",35)) oji<-c(rep("OJI",35)) area<-c(type_a,col,mar,oji) sabcolors <- c("MAR"="#8c510a", "OJI"="#dfc27d", "COL"="#35978f", "frugivore"="#e41a1c") # Junto vectores con dataset (cz_values) y transformo en numeric/factor y dataset. cz_values<-cbind(cz_values,type,module,area) cz_values$type<-as.factor(cz_values$type) cz_values$module<-as.factor(cz_values$module) cz_values$area<-as.factor(cz_values$area) cz_values$c<-as.numeric(cz_values$c) cz_values$z<-as.numeric(cz_values$z) cz_values<-as.data.frame(cz_values) ``` #Plotting -c -z values space ```{r} ggplot(data=cz_values, aes(x=c, y=z, label=row.names(cz_values), color =module )) +scale_color_manual(values=modcolors)+ geom_point() +geom_hline(yintercept=z.crit) + geom_vline(xintercept = c.crit) + theme_bw()+xlim(0, 0.9) +ylim(-1.3, 6) +geom_text_repel(aes(label=row.names(cz_values))) ggplot(data=cz_values, aes(x=c, y=z, label=row.names(cz_values), color =type, )) +scale_color_manual(values=type_color)+ geom_point() +geom_hline(yintercept=z.crit) + geom_vline(xintercept = c.crit) + theme_bw()+xlim(0, 0.9) +ylim(-1.3, 6) +geom_text_repel(aes(label=row.names(cz_values))) ggplot(data=cz_values, aes(x=c, y=z, label=row.names(cz_values), color =area )) +scale_color_manual(values=sabcolors)+ geom_point() +geom_hline(yintercept=z.crit) + geom_vline(xintercept = c.crit) + theme_bw()+xlim(0, 0.9) +ylim(-1.3, 6) +geom_text_repel(aes(label=row.names(cz_values))) ``` # ¿What occour with C153? ```{r} visits_mat<-read.csv("~/Documents/GitHub/Juniperus_mutualism_2022/data/visits_mat.csv",header=T,sep=";",dec=".",na.strings="NA") visits_mat<-column_to_rownames(visits_mat, var = "ind") visits_mat<-as.matrix(visits_mat) #heatmap(visits_mat) visits_mat #Is plant interacting with Genetta genetta. ``` Exploring the relationship between -c -z scores and fitness ```{r -c, -z prospection} #Parto dataset de cz_values para quedarme con las plantas solamente cz_values_plants<-(cz_values[13:117 ,1:2]) # Le añado a este dataset el producto de -c * -z como medida de rol topológico. cz_values_plants<-cbind(cz_values_plants,"cz_product"= (cz_values_plants$c * cz_values_plants$z)) hist(cz_values_plants$cz_product) plot(cz_values_plants$cz_product~dispersed_seeds_planttraits$dispersed_seeds) plot(cz_values_plants$c~dispersed_seeds_planttraits$dispersed_seeds) plot(cz_values_plants$z~dispersed_seeds_planttraits$dispersed_seeds) ``` ```{r} #dispersed_seeds_planttraits<-cbind(dispersed_seeds_planttraits,cz_values_plants) #dispersed_seeds_planttraits$c<-as.numeric(dispersed_seeds_planttraits$c) #dispersed_seeds_planttraits$z<-as.numeric(dispersed_seeds_planttraits$z) #mod_c<-lm(data =dispersed_seeds_planttraits,c~d1+h+cover+collapsed_cropsize+neighboiurhood_cropsize+n_neigh+long+dispersed_seeds) #summary(mod_c) #car::avPlots(mod_c,col.lines ="black") #jtools::plot_summs(mod_c) #mod_z<-lm(data =dispersed_seeds_planttraits,z~d1+h+cover+collapsed_cropsize+neighboiurhood_cropsize+n_neigh+long+dispersed_seeds) #jtools::plot_summs(mod_z) ##summary(mod_z) #car::avPlots(mod_z,col.lines ="black") #relaimpo::calc.relimp(mod_z,type = "lmg") ``` Splitting networks by module to analyse within-module structure ```{r} to_cut<-cbind(visits_mat,"module" =plants_105$module) to_cut<-as.data.frame(to_cut) to_cut$module<-as.factor(to_cut$module) require(dplyr) mod_A<-filter(to_cut, to_cut$module == "A") mod_B<-filter(to_cut, to_cut$module == "B") mod_C<-filter(to_cut, to_cut$module == "C") mod_A<-as.matrix(mod_A[,c(1,11)]) mod_B<-mod_B[,c(2,3,4)] mod_C<-mod_C[,c(5,6,7,8,9,10,12)] mod_A<-as.data.frame(mod_A) mod_B<-as.data.frame(mod_B) mod_C<-as.data.frame(mod_C) #A mod_A$tur_phi<-as.numeric(mod_A$tur_phi) mod_A$gen_gen<-as.numeric(mod_A$gen_gen) #B mod_B$vul_vul<-as.numeric(mod_B$vul_vul) mod_B$tur_mer<-as.numeric(mod_B$tur_mer) mod_B$eri_rub<-as.numeric(mod_B$eri_rub) #C mod_C$tur_ili<-as.numeric(mod_C$tur_ili) mod_C$syl_mel<-as.numeric(mod_C$syl_mel) mod_C$syl_atr<-as.numeric(mod_C$syl_atr) mod_C$ory_cun<-as.numeric(mod_C$ory_cun) mod_C$cya_cya<-as.numeric(mod_C$cya_cya) mod_C$tur_tor<-as.numeric(mod_C$tur_tor) mod_C$mel_mel<-as.numeric(mod_C$mel_mel) plotweb(mod_A, method = "normal") plotweb(mod_B, method = "normal") plotweb(mod_C, method = "normal") visweb(mod_A,type = "nested") visweb(mod_B,type = "nested") visweb(mod_C,type = "nested") networklevel(mod_A, index="weighted NODF", level="both", weighted=TRUE) networklevel(mod_B, index="weighted NODF", level="both", weighted=TRUE) networklevel(mod_C, index="weighted NODF", level="both", weighted=TRUE) ``` General centrality in relationship with dispersed seeds ```{r} #Estimating specieslevel-based metrics species_metrics_animals<-specieslevel(visits_mat, index=c("betweenness","degree", "species strength", "closeness"), level="higher") species_metrics_plants<-specieslevel(visits_mat, index=c("betweenness","degree", "species strength", "closeness"), level="lower") hist(species_metrics_plants$degree) hist(species_metrics_plants$species.strength) plot(species_metrics_plants$weighted.betweenness) plot(species_metrics_plants$weighted.closeness) #row.names(species_metrics_plants) #row.names(dispersed_seeds_planttraits) dispersed_seeds_planttraits<-cbind(dispersed_seeds_planttraits,species_metrics_plants) dispersed_seeds_planttraits<-as.data.frame(dispersed_seeds_planttraits) # Fitting model mod_s<-lm(data =dispersed_seeds_planttraits,weighted.closeness~dispersed_seeds) summary(mod_s) residuals(mod_s) plot(mod_s) # PLOTTING names<-row.names(dispersed_seeds_planttraits) p2 <- ggplot(dispersed_seeds_planttraits, aes(x=weighted.closeness, y=dispersed_seeds, label = names)) p2<-p2 +geom_text_repel() p2<-p2 + geom_point(aes(color = module, size = 0.5))+ scale_color_manual(values = c("A"="#b2182b", "B"="#f4a582", "C"="#bababa")) p2 + stat_smooth(method = "lm", col = "grey") +theme_light() correlation <- cor.test(dispersed_seeds_planttraits$weighted.closeness, dispersed_seeds_planttraits$dispersed_seeds, method = 'pearson') correlation ```