###Code under R environment

###################
# Code for: "Delimiting the boundaries of sesamoid identities under the network theory framework"
# (Fontanarrosa, Fratani & Vera, 2020)  

library(igraph)
library(onewaytests)
library(brainGraph)

SESRAW <- read.csv(choose.files(), header = F, sep=";") # Read the attached file Supplemental Table S4 ("TableS4_networkdata.csv")

AM<-get.adjacency(graph.edgelist(as.matrix(SESRAW), directed=FALSE))
gSES <- graph_from_adjacency_matrix(AM, mode = "undirected", weighted = TRUE, diag = FALSE)
E(gSES)$width <- E(gSES)$weight 

# Network parameters
V(gSES)	# Node number
E(gSES)	# Edge (links) number
edge_density(gSES, loops = FALSE)	# Density
average.path.length(gSES, directed=TRUE, unconnected=TRUE)	# Average shortest paths 
transitivity(gSES, type = "undirected", vids=V(graph))	# Global clustering coefficient

# Centrality indicators
grado <- degree(gSES)	# Degree
clos <- closeness(gSES, vids = V(gSES), normalized=TRUE)	# Closeness
inter <- betweenness(gSES)	# Betweenness
eigen <- eigen_centrality(gSES)$vector # Eigen_centrality

# Shortest Paths
CMC <- shortest.paths(gSES, V(gSES), V(gSES))

# Transitivity of each node
trans <- transitivity(gSES, type = "local", vids=V(gSES))           


tkplot(gSES) # Interactive window for network visualization.

##### Centrality values comparisons among categories
# Degree
grado<-read.csv(choose.files(), header=T, sep=";") # Choose the attached file Supplemental Table S5  ("TableS5_degree.csv"). 

# mean calculations by categories
mean((subset(grado,grado$Cat=="SE"))$deg)
mean((subset(grado,grado$Cat=="SG"))$deg)
mean((subset(grado,grado$Cat=="NS"))$deg)

# Kruskal Wallis test comparing among all categories 

kruskal.test(deg~Cat, data=grado)
# Post-hoc Wilcoxon pairwise comparisons among all categories

pairwise.wilcox.test(grado$deg, grado$cat, p.adjust.method="BH")

# Boxplot - degree
boxplot(grado$deg~grado$Cat, ylab="Degree", xlab="Node skeletal category", col=c("green", "blue", "red"),  boxwex=.21)

# Closeness
Clos<-read.csv(choose.files(), header=T, sep=";")  # Choose the attached file Supplemental Table S6 ("TableS6_closeness.csv").
# mean calculations by categories
mean((subset(Clos,Clos$Cat=="SE"))$clos)
mean((subset(Clos,Clos$Cat=="SG"))$clos)
mean((subset(Clos,Clos$Cat=="NS"))$clos)

# Kruskal Wallis  test comparing among all categories 
kruskal.test(clos~Cat, data=Clos)

# Post-hoc Wilcoxon pairwise comparisons among all categories
pairwise.wilcox.test(Clos$clos, Clos$Cat, p.adjust.method="BH")

# Boxplot - closeness
boxplot(Clos$clos~Clos$Cat, ylab="Closeness", xlab="Node skeletal category", col=c("green", "blue", "red"),  boxwex=.21)

# Betweenness
Bet<-read.csv(choose.files(), header=T, sep=";")# Read the file Supplemental Table S7 ("TableS7_betweenness.csv").
# mean calculations by categories
mean((subset(Bet,Bet$Cat=="SE"))$bet)
mean((subset(Bet,Bet$Cat=="SG"))$bet)
mean((subset(Bet,Bet$Cat=="NS"))$bet)
# Kruskal Wallis test comparing among all categories 
kruskal.test(bet~Cat, data=Bet)
#Post-hoc Wilcoxon pairwise comparisons among all categories
pairwise.wilcox.test(Bet$bet, Bet$Cat, p.adjust.method="BH")
# Boxplot - Betweenness
boxplot(Bet$bet~Bet$Cat, ylab="Betweenness", xlab="Node skeletal category", col=c("green", "blue", "red"))

# Eigen-centrality
Eig<-read.csv(file=choose.files(), header=T, sep=";") # Read the attached file Supplemental Table S8 ("TableS8_eigen_centrality.csv")
# mean calculations by categories
mean((subset(Eig,Eig$Cat=="SE"))$eig)
mean((subset(Eig,Eig$Cat=="SG"))$eig)
mean((subset(Eig,Eig$Cat=="NS"))$eig)

# Kruskal Wallis test comparing among all categories
kruskal.test(eig~Cat, data=Eig)
# Post-hoc Wilcoxon pairwise comparisons among all categories
pairwise.wilcox.test(Eig$eig, Eig$Cat, p.adjust.method="BH")
# Boxplot - Eigen-centrality
boxplot(Eig$eig~Eig$Cat, ylab="Eigen-centrality", xlab="Node skeletal category", col=c("green", "blue", "red"))

# Base for Figure 1 -----------------
modules <- read.csv(file=choose.files(), header=T, sep=";")  # Read the attached file Supplemental Table S9 ("TableS9_modules.csv")
head(modules) # Modules and colors are assigned 

## Plot
gSES <- graph.edgelist(as.matrix(SESRAW), directed = F)
V(gSES)$color <- modules$color[match(V(gSES)$name,modules$Element)]

plot.igraph(gSES, vertex.label = NA, vertex.size=2.5, layout=layout_with_graphopt)
tkplot(gSES) # Interactive window for network visualization.