ISMEJC_data_assembly.Rmd

---
title: "R Notebook"
output: html_notebook
---


# Loading the data

Data from the analysis of mycelial pictures using the network app from Mark Fricker
```{r}

rm(list = ls())
library(readxl)
library(tidyverse)
data_location<-"ISMEJ_C_data//"#This is the latest data as it comes from the most 

#Locating the name of the new datafile
temp = list.files(path=data_location,pattern="*-results.xlsx")
temp<-paste(data_location,temp,sep = "")


#Edges
Edge_Traits = lapply(temp,function(x){read_excel(x,sheet = "Edges")} )
                

#Nodes
Node_Traits = lapply(temp,function(x){read_excel(x,sheet = "Nodes")} )

#Summary_Tables
summary_Traits = lapply(temp,function(x){read_excel(x,sheet = "Summary table")} )


summary_Traits_table_AZB <- bind_rows(summary_Traits)
summary_Traits_table_AZB$name_col <- temp
summary_Traits_table_AZB$name_col <- gsub(data_location,"",summary_Traits_table_AZB$name_col)
summary_Traits_table_AZB$name_col <- gsub("\\-results\\.xlsx","",summary_Traits_table_AZB$name_col)

#Ordered robustness
Ordered_robustness <- lapply(temp,function(x){read_excel(x,sheet = "Ordered Robustness")})

#Spatial robustness
Spatial_robustness <- lapply(temp,function(x){read_excel(x,sheet = "Spatial Robustness")})

#Random robustness
Random_robustness <- lapply(temp,function(x){read_excel(x,sheet = "Random Robustness")})

```

# Data preparation

## Correcting values

The value of resistance (Resistance_2ave) for microscopic fungi has to be updated. This is because, as we mentioned in the materials and methods of the paper, resistance for microscopic fungi should follow: 16*(Length/Width^4).

The current value of resistance is based on the relationship for macroscopic fungi which follows: 4*(Length/Width^4).  The relationship for microscopic fungi is already calculated in Resistance_4ave. Thus for simplicity in the analysis, I will just replaced all the values of Resistance_2ave with Resistance_4ave ONLY for the microfungi. In that way, the rest of the code follows unaltered.

```{r}
#Colonies 1 to 15 correspond to microfungi (ascos and zygos). One can see that with commented out code: sapply(Edge_Traits[1:26],function(x){unique(x$name)}) Colony 26 is a basidio alerady

Edge_Traits[1:25]<-
lapply(
Edge_Traits[1:25],function(x){x$Resistance_2ave<-x$Resistance_4ave
x
})
```

Calculating branching angle

The branching angle correspond the  smallest angle calcuated by the FungalNetworks App out of these three: "node_Omin_Omaj","node_Omid_Omaj","node_Omin_Omid"


Note: This code has to be run in the console only!!! It turns out rmarkdown gets caught up in a loop
```{r, eval=FALSE}

Node_Traits_c<-
lapply(Node_Traits,function(prueba){
  
  prueba$min_angle<-apply(
        prueba[,c("node_Omin_Omaj","node_Omid_Omaj","node_Omin_Omid")],1,
        function(x){y<-min(x,na.rm=T)})

  prueba$min_angle[which(prueba$min_angle=="Inf")]<-NA
  prueba
  })

all(Node_Traits_c[[1]]$node_Strength==Node_Traits[[1]]$node_Strength)
Node_Traits<-Node_Traits_c;rm(Node_Traits_c)
```


## Transforming it into igraph object

```{r, echo=FALSE, message=FALSE, warning=FALSE}
#Transforming it into igraph object
library(igraph)
library(tidyverse)

colonies_ntwk<-

mapply(function(x,z){
    y<-graph_from_edgelist(as.matrix(x[,c("EndNodes_1","EndNodes_2")]),directed = F)
    E(y)$name<-x$name
    E(y)$weight<-x$Resistance_2ave
    E(y)$length<-x$Length
    E(y)$width<-x$Width
    E(y)$Volume <- x$Volume
    
    V(y)$Degrees<-degree(y)
    V(y)$Accessibility<-distances(y,v=as.numeric(V(y)[Degrees==max(V(y)$Degrees)])
    
    
                                                                )
    V(y)$node_ID<-z$node_ID
    #V(y)$angle<-z$node_Omin_Omid
    V(y)$angle<-z$min_angle
    
    E(y)$type<-"Main"
    
    E(y)$e_distance<-(x$Length/x$Tortuosity)
    
    E(y)[incident(y,
              as.numeric(V(y)[Degrees==max(V(y)$Degrees)])
                    )]$type<-"Inoculum"
    
    y
    },Edge_Traits,Node_Traits,SIMPLIFY = F)

#saveRDS(colonies_ntwk,"colonies_ntwk.RDS")
```


Saving coordinates as a separate file
```{r, echo=FALSE, message=FALSE, warning=FALSE}
#Spatial location
spatial.data<-lapply(Node_Traits,function(x){x[,c("node_ID","node_X_pix","node_Y_pix","name")]})
    
spatial.data<-lapply(spatial.data,function(l){as.matrix(l[,c(2,3)])})
spatial.data<-lapply(spatial.data,function(l){norm_coords(l, ymin=-1, ymax=1, xmin=-1, xmax=1)})

#saveRDS(spatial.data,"spatial.data.RDS")
```

Calculating the summary (mean) values for each colony (hyphae traits like length and width)

```{r, echo=FALSE, message=FALSE, warning=FALSE}
#Calculating the summary values for each colony

colony_sum<-do.call("rbind",

lapply(
colonies_ntwk,function(x){
  
  w<-data.frame(
          name_col=unique(E(x)$name),
          alpha_coeff=(ecount(x)-vcount(x)+1)/(2*vcount(x)-5),
          beta_coeff=ecount(x)/vcount(x),
          gamma_coeff=ecount(x)/(3*vcount(x)-6),
          Root_eff=mean(1/V(x)[Accessibility>0]$Accessibility),
          Reff_tip=mean(1/V(x)[Degrees==1&Accessibility>0]$Accessibility))
  
    y<-subgraph.edges(x,E(x)[type=="Main"])
    
    E(y)$hyphae<-"main"
    E(y)[inc(V(y)[Degrees==1])]$hyphae<-"tip"
    
    z<-data.frame(
        name_col=unique(E(y)$name),
        Hyphal_length=mean(log10(E(y)$length)),
        Hyphal_number=length(E(y)$name),
        Hyphal_tip_width=mean(log10(E(y)[hyphae=="tip"&width>0]$width)),
        Hyphal_main_width=mean(E(y)[hyphae=="main"&width>0]$width),#The unstransformed looks better
        Hyphal_tip_access=mean(V(y)[Degrees==1]$Accessibility),
        Hyphal_angle=mean(V(y)[which(angle>0)]$angle),
        
        Mycelia_length=sum(E(y)$length),
        Mycelia_Volume=sum(E(y)$Volume)
        
            )
    
    zw<-left_join(z,w)
    zw
    
    })
)

```


Calculating global efficiency

```{r}
library(parallel)
library(doParallel)

global_eff<-
function(g){

get_eff <- function(i){return((1/(length(V(g)) - 1))*sum(1/distances(g, V(g)[i])[-i]))}
no_cores <- detectCores() - 1 
cl       <- makeCluster(no_cores)
registerDoParallel(cl)  
result <- foreach(i = seq_along(V(g)), .combine = c, .packages = "igraph") %dopar% get_eff(i)
stopCluster(cl)
rm(cl)
Geff<-mean(result);Geff
}

Real_summary_Geff1<-#it took 1 hour in my laptop (16GB) and 25 min in the FU desktop (64 GB)
  sapply(colonies_ntwk,global_eff)

saveRDS(Real_summary_Geff1,"ISMEJ_C_output//Real_summary_Geff.RDS")
  
```


```{r}
Real_summary_Geff<-readRDS("ISMEJ_C_output//Real_summary_Geff.RDS")
colony_sum$Geff<-Real_summary_Geff
```

Calculating (weighted) Minimum spanning trees based on the dimesnions of each of the networks. The weights are given by resistance_2

```{r,echo=FALSE, message=FALSE, warning=FALSE}
#Weighted MST (by accessibility)
Access_mst<-
        lapply(colonies_ntwk,
               mst)

#saveRDS(Access_mst,"Access_mst.RDS")
```


Calculating summary values (hypal traits, coefficients and efficiencies ) for the minimum spanning trees (MST´s)

```{r,echo=FALSE, message=FALSE, warning=FALSE}
#Summarizng the MST´s unweighted_mst, edistance_mst, Access_mst_sum

Access_mst_sum<-do.call("rbind",
    lapply(
    Access_mst,function(x){
      
      w<-data.frame(
          name_col=unique(E(x)$name),
          alpha_coeff=(ecount(x)-vcount(x)+1)/(2*vcount(x)-5),
          beta_coeff=ecount(x)/vcount(x),
          gamma_coeff=ecount(x)/(3*vcount(x)-6),
          Root_eff=mean(1/V(x)[Accessibility>0]$Accessibility),
          Reff_tip=mean(1/V(x)[Degrees==1&Accessibility>0]$Accessibility))
      
        y<-subgraph.edges(x,E(x)[type=="Main"])
        
        E(y)$hyphae<-"main"
        E(y)[inc(V(y)[Degrees==1])]$hyphae<-"tip"
        
        z<-data.frame(
            name_col=unique(E(y)$name),
            Hyphal_length=mean(E(y)$length),
            Hyphal_number=length(E(y)$name),
            Hyphal_tip_width=mean(E(y)[hyphae=="tip"]$width),
            Hyphal_main_width=mean(E(y)[hyphae=="main"]$width),
            Hyphal_tip_access=mean(V(y)[Degrees==1]$Accessibility),
            Hyphal_angle=mean(V(y)[which(angle>0)]$angle),
            
            Mycelia_length=sum(E(y)$length),
            Mycelia_Volume=sum(E(y)$Volume)#,
            
            #Now the indexes and efficiencies
            # alpha_coeff=(ecount(x)-vcount(x)+1)/(2*vcount(x)-5),
            # beta_coeff=ecount(x)/vcount(x),
            # gamma_coeff=ecount(x)/(3*vcount(x)-6),
            # Root_eff=mean(1/V(x)[Accessibility>0]$Accessibility),
            # Reff_tip=mean(1/V(x)[Degrees==1&Accessibility>0]$Accessibility)
        )
        z<-left_join(z,w)
        
        })
    )
```


For the global efficiency, It was more computing demanding. I save the results for not having to do it again

```{r}
MST_summary_Geff1<-#do.call(rbind,#took 1 hour too in my laptop
sapply(Access_mst,global_eff)
  
saveRDS(MST_summary_Geff1,"ISMEJ_C_output//MST_summary_Geff.RDS")
```

```{r}
MST_summary_Geff<-readRDS("ISMEJ_C_output//MST_summary_Geff.RDS")
Access_mst_sum$Geff<-MST_summary_Geff
```

### Biding the real and MST networks

```{r, echo=FALSE, message=FALSE, warning=FALSE}
colony_sum$Network<-"Real"
Access_mst_sum$Network<-"Resistance_MST"

all_data<-bind_rows(colony_sum,Access_mst_sum)

#Adding the species names
all_data$Species<-NA
all_data$Species[grep("C34",all_data$name_col)]<-"Mortierella elongata"
all_data$Species[grep("C35",all_data$name_col)]<-"Umbelopsis isabellina"
all_data$Species[grep("DF19",all_data$name_col)]<-"Mortierella alpina"
all_data$Species[grep("DF25",all_data$name_col)]<-"Mortierella elongata2"
all_data$Species[grep("DF56",all_data$name_col)]<-"Mucor fragilis"
all_data$Species[grep("M",all_data$name_col)]<-"Mortierella alpina2"

all_data$Species[grep("DF9",all_data$name_col)]<-"Alternaria sp"
all_data$Species[grep("C41",all_data$name_col)]<-"Fusarium redolens"
all_data$Species[grep("FOX",all_data$name_col)]<-"Fusarium oxysporum"
all_data$Species[grep("DF32",all_data$name_col)]<-"Fusarium solani"

all_data$Species[grep("Pi",all_data$name_col)]<-"Phallus impudicus"
all_data$Species[grep("Pv",all_data$name_col)]<-"Phanerachaete ventulina"
all_data$Species[grep("Rb",all_data$name_col)]<-"Resinicium bicolor"
all_data$Species[grep("Hf",all_data$name_col)]<-"Hypholoma fasiculare"

#Adding area data
all_data<-left_join(all_data,
summary_Traits_table_AZB[,c("name_col","summary_mycelial_area")])


#Then I can just calculated edge density for all (including the MSt´s)

all_data$Hyphal_density<-all_data$Hyphal_number/all_data$summary_mycelial_area

#Changin the name of the area so it make more sense
names(all_data)[which(names(all_data)=="summary_mycelial_area")]<-"Mycelial_area"

```

### Scaling values from real networks by MST values

```{r}
l<-length(all_data$name_col)


all_data_scaled<-all_data[c(1:(l/2)),]

all_data_scaled$Mycelia_length_MST<-all_data[c(1:(l/2)),]$Mycelia_length/all_data[c(((l/2)+1):l),]$Mycelia_length

all_data_scaled$Mycelia_Vol_MST<-all_data[c(1:(l/2)),]$Mycelia_Volume/all_data[c(((l/2)+1):l),]$Mycelia_Volume

#all_data_scaled$alpha_coeff_scaled<-all_data[c(1:32),]$alpha_coeff/all_data[c(33:64),]$alpha_coeff
#all_data_scaled$beta_coeff_scaled<-all_data[c(1:32),]$beta_coeff/all_data[c(33:64),]$beta_coeff
all_data_scaled$Root_eff_l_scaled<-all_data[c(1:(l/2)),]$Root_eff/sqrt(all_data[c(1:(l/2)),]$Mycelial_area)
all_data_scaled$Reff_tip_l_scaled<-all_data[c(1:(l/2)),]$Reff_tip/sqrt(all_data[c(1:(l/2)),]$Mycelial_area)
all_data_scaled$Geff_MST<-all_data[c(1:(l/2)),]$Geff/all_data[c(((l/2)+1):l),]$Geff

```

Adding phylum data
```{r}
all_data_scaled$phylum<-NA
all_data_scaled$phylum[grep("DF9|C41|FOX|DF32",all_data_scaled$name_col)]<-"Ascomycota"
all_data_scaled$phylum[grep("C34|C35|DF19|DF25|DF56|M",all_data_scaled$name_col)]<-"Zygomycetous"
all_data_scaled$phylum[grep("Pi|Pv|Rb|Hf",all_data_scaled$name_col)]<-"Basidiomycota"

#saveRDS(all_data_scaled,"all_data_scaled.RDS")
```


# Adding robustness measures (For ascos and zygos)

## Ordered robustness
```{r}
#adding a new column name
Ordered_robustness<-
lapply(
Ordered_robustness,function(x){
  x$dir_method<-tolower(paste(x$direction,x$method,sep = "_"))
  x
})


# 1. Fifty percent cut-off
Ordered_50<-do.call("rbind",
                   sapply(Ordered_robustness,function(datos){
                     lapply(split(datos,datos$dir_method),
                            function(x){y<-data.frame(name=unique(x$name),
                                                      dir_method=unique(x$dir_method),
                                                      Fifty_mark= #approx(x$robustness,x$removed_percent,xout=50)$y)
approx(x$removed_percent,x$robustness)$x[
which.min(abs(approx(x$removed_percent,x$robustness)$y-50))])                                                        
                            })
                   })
)

Ordered_50<-
Ordered_50 %>% 
pivot_wider(values_from = Fifty_mark,names_from=dir_method)

names(Ordered_50)[-1]<-paste(names(Ordered_50)[-1],"robustness",sep = "_")
```


## Spatial robustness
```{r}

names(Spatial_robustness)<-sapply(Spatial_robustness,function(x){unique(x$name)})

# 1. Fifty percent cut-off

Spatial_50<-do.call("rbind",
                    #sapply(Random_robustness,function(datos){
                    lapply(Spatial_robustness,function(x){
                      #x<-aggregate(robustness~name*removed_percent,datos,mean)
                      z<-data.frame(name=unique(x$name),
                                    Fifty_mark= #approx(x$robustness,x$removed_percent,xout=50)$y)
approx(x$removed_percent,x$robustness)$x[
which.min(abs(approx(x$removed_percent,x$robustness)$y-50))])                      
                      z}
                    ))

Spatial_50$Fifty_mark[which(is.na(Spatial_50$Fifty_mark))]<-99
rownames(Spatial_50)<-NULL
names(Spatial_50)[2]<-"spatial_robustness"
```


## Random robustness
```{r}

names(Random_robustness)<-sapply(Random_robustness,function(x){unique(x$name)})
# 1- Fifty percent mark
Random_50<-do.call("rbind",
                   sapply(Random_robustness,function(datos){
                     lapply(split(datos,as.factor(datos$rep)),
                            function(x){y<-data.frame(name=unique(x$name),
                                                      rep=unique(x$rep),
                                                      Fifty_mark= #approx(x$robustness,x$removed_percent,xout=50)$y)
approx(x$removed_percent,x$robustness)$x[
which.min(abs(approx(x$removed_percent,x$robustness)$y-50))])                            
                            })
                   })
);Random_50<-do.call("rbind",Random_50)

Random_50 <-aggregate(Fifty_mark~name,Random_50,mean)
names(Random_50)[2]<-"random_robustness"
```

# Merging them all
```{r}
#Just standardizing the name column because it is what will be used to left_join
names(Ordered_50)[1]<-"name_col"
names(Random_50)[1]<-"name_col"
names(Spatial_50)[1]<-"name_col"

#Now merging them
robustness<-
left_join(Ordered_50,Random_50,by="name_col")

robustness<-
left_join(robustness,Spatial_50,by="name_col")

robustness$asc_width_robustness[which(is.na(robustness$asc_width_robustness))]<-99

#saveRDS(robustness,"robustness.RDS")
```

Saving data for analysis
```{r}
saveRDS(colonies_ntwk,"ISMEJ_C_output//colonies_ntwk.RDS")
saveRDS(Access_mst,"ISMEJ_C_output//Access_mst.RDS")
saveRDS(all_data_scaled,"ISMEJ_C_output//all_data_scaled.RDS")
saveRDS(all_data,"ISMEJ_C_output//all_data.RDS")
saveRDS(robustness,"ISMEJ_C_output//robustness.RDS")
saveRDS(spatial.data,"ISMEJ_C_output//spatial.data.RDS")

saveRDS(Edge_Traits,"ISMEJ_C_output//Edge_Traits.RDS")
saveRDS(Node_Traits,"ISMEJ_C_output//Node_Traits.RDS")


#summary values from Mark
saveRDS(summary_Traits_table_AZB,"ISMEJ_C_output//summaryTraits_Mark_AZB.RDS")

```