#################################################################################################################################
###  Multiannual co-occurrences between Saimaa ringed seal individuals suggest social associations during the moulting season ###
#################################################################################################################################
									
#####################################
### Network Construction (Fig.2) ####
#####################################

### Load packages ###
library(png)
library(igraph)

### Individual Info###
ID <- read.table(file.choose(), header=T)
#File: "Individual_Info_All_Basins.txt"

### Factor correction ###
ID$ID <- as.character(ID$ID)
ID$Sex <- as.character(ID$Sex)
ID$Basin <- as.character(ID$Basin)

### Matrix ###
matrix <- read.table(file.choose(), header=TRUE, row.names=1, check.names=FALSE)
m <- as.matrix(matrix)
diag(m)=0
#File: "Matrix_Cooccurences.txt"
#The matrix was build with SOCPROG2.10 with the parameters mention in the Methods section of the paper.

### Keep the edge value ###
net <- graph.adjacency(m, mode="undirected", weighted=TRUE)

### Network Overview ###
plot(net)

#Caluclate individuals' degree
ID$Degree <- rowSums(m !=0)

### If you want the node to have different shapes ###
### By default, igraph has only circle and square. Next function to create triangle shapes###
### triangle vertex shape ###
mytriangle <- function(coords, v=NULL, params) {
  vertex.color <- params("vertex", "color")
  if (length(vertex.color) != 1 && !is.null(v)) {
    vertex.color <- vertex.color[v]
  }
  vertex.size <- 1/200 * params("vertex", "size")
  if (length(vertex.size) != 1 && !is.null(v)) {
    vertex.size <- vertex.size[v]
  }

  symbols(x=coords[,1], y=coords[,2], bg=vertex.color,
          stars=cbind(vertex.size, vertex.size, vertex.size),
          add=TRUE, inches=FALSE)
}
add_shape("triangle", clip=shape_noclip,
                 plot=mytriangle)
				 
				 
### Give shape to Sex ###
ID$Shape<- "grey"
ID$Shape[which(ID$Sex == "M")] <- "square"
ID$Shape[which(ID$Sex == "F")] <- "circle"
ID$Shape[which(ID$Sex == "U")] <- "triangle"

#Give symbol Colour according to lake basin
ID$Colour <- "grey"
ID$Colour[which(ID$Basin == "PS")] <- "#9e0142" #bordeaux
ID$Colour[which(ID$Basin == "PEV")] <- "#d53e4f" #brick
ID$Colour[which(ID$Basin == "KV")] <- "#f46d43" #dark orange
ID$Colour[which(ID$Basin == "JV")] <- "#fdae61" #orange
ID$Colour[which(ID$Basin == "HV")] <- "#fee08b" #mustard
ID$Colour[which(ID$Basin == "PV")] <- "#abdda4" #green
ID$Colour[which(ID$Basin == "PUV")] <- "#e6f598" #light green
ID$Colour[which(ID$Basin == "KS")] <- "#66c2a5" #aqua
ID$Colour[which(ID$Basin == "LL")] <- "#3288bd" #blue
ID$Colour[which(ID$Basin == "ES")] <- "#5e4fa2" #purple


#Edge label
Elabel <- E(net)$weight

#Check the network with the edge values
plot(net, vertex.label=NA, vertex.color=ID$Colour, vertex.label=NA, vertex.size=ID$Degree*3, vertex.shape=ID$Shape,edge.color="black", edge.width=E(net)$weight, edge.label=Elabel, edge.label.cex = 3)

### Adjust the layout with tkplot ###
### Change the vertex position manualy and save coordinates ###
tkplot(net, vertex.label=NA, vertex.color=ID$Colour, vertex.label=NA, vertex.size=ID$Degree*3, vertex.shape=ID$Shape,edge.color="black", edge.width=E(net)$weight, edge.label=Elabel, edge.label.cex = 3)

### Tabulate the coordinates ###
tk_coords(id, norm=FALSE) -> Plot_coord
### Export coordinates ###
write.table(Plot_coord, file="Plot_coord.txt", sep="\t")

### Import the coordinates ###
Plot_coord <- read.table(file.choose(), header=F)
Plot_coord <- as.matrix(Plot_coord)
#File: "Plot_coord.txt"
#These are the coordinates used for Fig. 2 of the paper

###Final Plot###
plot(net, vertex.label=NA, vertex.color=ID$Colour, vertex.label=NA, vertex.size=ID$Degree*3, vertex.shape=ID$Shape,edge.color="black", edge.width=E(net)$weight, edge.label.cex = 3, layout=Plot_coord)

###With edge values
plot(net, vertex.label=NA, vertex.color=ID$Colour, vertex.label=NA, vertex.size=ID$Degree*3, vertex.shape=ID$Shape,edge.color="black", edge.width=E(net)$weight,edge.label=Elabel, edge.label.cex = 3, layout=Plot_coord)


### Export Final Plot ###
tiff("Network_All_Basins.tiff", width=12, height=12, units='in', res=900)
plot(net, vertex.label=NA, vertex.color=ID$Colour, vertex.label=NA, vertex.size=ID$Degree*3, vertex.shape=ID$Shape,edge.color="black", edge.width=E(net)$weight, edge.label.cex = 3, layout=Plot_coord)
dev.off()