#To generate bipartite networks and calculate networks metrics using the bipartite package. Do the same in "SAVANNA" and "TOTAL" files.
install.packages("bipartite")
library(bipartite)
rainforestnet<-read.delim("RAINFOREST.txt", row.names = 1)
networklevel(rainforestnet)
specieslevel(rainforestnet) 
plotweb(rainforestnet, col.high="darkgoldenrod1", col.low="chartreuse4", arrow= "down", text.rot = 90) 

###Add circular graph using igraph package. Do the same in "RAINFOREST2" and "SAVANNA2". Thanks to Spiesman Lab (https://spiesmanecology.com/2017/04/30/network-vis1/) for the R script.
install.packages("igraph")
library(igraph)
net<-read.delim("TOTAL2.txt", row.names = 1)
i_net = graph_from_incidence_matrix(net, weight=T)
plot(i_net)
plot(i_net, layout=layout.circle)
V(i_net)$frame.color = "white"
plot(i_net,layout=layout.circle)

#For degree centrality:
deg = centr_degree(i_net, mode="all")
V(i_net)$size = 5*sqrt(deg$res)
plot(i_net, layout=layout.circle)
polcol = rep("gray50",dim(net)[2])
plantcol = rep("black",dim(net)[1])
clrs = rbind(as.matrix(plantcol),as.matrix(polcol))
V(i_net)$color = clrs
plot(i_net,layout=layout.circle)
E(i_net)$width = E(i_net)$weight/15
plot(i_net,layout=layout.circle)

#Identify modules for plants and pollinators
i_net.com = fastgreedy.community(i_net)
plantmod = i_net.com$membership[1:dim(net)[1]]
pollmod = i_net.com$membership[(dim(net)[1]+1):(length(i_net.com$membership))]

#Re-order the network matrix rows and columns by module ID
onet = net[order(plantmod),] 
onet = t(t(onet)[order(pollmod),])
i_net = graph_from_incidence_matrix(onet, weight=T)

#Add a color palette
install.packages("wesanderson")
library(wesanderson)
com.mem = rbind(as.matrix(sort(plantmod)),as.matrix(sort(pollmod)))
colrs = wes_palette(name="Darjeeling1", 
                    n=max(i_net.com$mem),
                    type="continuous")

#Specify vertex and edge parameters
deg = centr_degree(i_net, mode="all")
V(i_net)$size=5*sqrt(deg$res)
V(i_net)$size=20
V(i_net)$frame.color="white"
V(i_net)$color=colrs[com.mem]
V(i_net)$label.color="black"
V(i_net)$label.font=2
V(i_net)$label.cex=1.5
E(i_net)$width=E(i_net)$weight/3
plot(i_net, layout=layout_in_circle, edge.color="grey15")

#To generate null matrices (r2dtable), calculate z scores and z-test. Do the same in H2', NODF and WNODF in "SAVANNA2" and "TOTAL2" files.
rainforestnet <-read.delim("RAINFOREST2.txt", row.names = 1)
networklevel(rainforestnet) 
obs_c <- unlist(networklevel(rainforestnet, index="connectance"))
obs_c
null1<-nullmodel(rainforestnet, N=1000, method= 1)
null <- unlist(sapply(null1, networklevel, index="connectance")) 
plot(density(null), xlim=c(min(obs_c, min(null)), max(obs_c, max(null))), 
     main="comparison of observed with null model Patefield")
abline(v=obs_c, col="red", lwd=2)   
praw <- sum(null>obs_c) / length(null)
praw
ifelse(praw > 0.5, 1-praw, praw)    # P-value
mean_c<-mean(null)
sd_c_rd2<-sd(null)
z<-(obs_c-mean_c)/sd_c_rd2
z
p <- 2*pnorm(-abs(z))
p

#To calculate network modularity using QuanBiMo. Do the same in "SAVANNA" and "TOTAL" files.
install.packages("bipartite")
library(bipartite)
rainforestnet <-read.delim("RAINFOREST.txt", row.names = 1)
modulesrainforest=computeModules(web= rainforestnet, steps=1e8)
printoutModuleInformation(modulesrainforest)
moduloslista=listModuleInformation(modulesrainforest)
plotModuleWeb(modulesrainforest)

#To calculate null models in network modularity, compute z scores and z-test. Do the same in "SAVANNA" and "TOTAL" files.
nulls <- nullmodel(rainforestnet, N=100, method=3) 
modules.nulls <- sapply(nulls, computeModules)
like.nulls <- sapply(modules.nulls, function(x) x@likelihood)
z <- (0.2920 - mean(like.nulls))/sd(like.nulls)
p <- 2*pnorm(-abs(z))

#To calculate phenotypic floral integration using PHENIX Package. Do the same for all plant species.
install.packages("PHENIX")
library(PHENIX)
aechbract<-read.delim("AECHBRACT_INT.txt")
#Partial correlations
cor.par(aechbract[,1:3], aechbract[,4], trait.names=TRUE) 
traits<-na.exclude(aechbract) # NA manually removed
cor.par(traits[,1:3], traits[,4], trait.names=TRUE)
#Phenotipic integration index
pint (aechbract[,1:3])