---
title: "First_workflow"
author: "Bingrui Chen"
date: "2024-06-01"
output:
  html_document:
    toc: true
    toc_depth: 4
    toc_float: true
    anchor_sections: true
    code_download: true
    highlight: textmate
    theme: journal
  pdf_document:
    toc: true
    toc_depth: '4'
---

```{r setup, include=FALSE}
#Set working path
knitr::opts_knit$set(
  root.dir = "#:/#/#/#/#"#pathway of your files
)

```

#1 R package loading

```{r}

library(dplyr)
library(readxl)
library(ggplot2)
library(vegan)
library(ggrepel)
library(microeco)
library(rcompanion)
library(gridExtra)
library(Biostrings)
library(readr)
library(tibble)
library(tidyverse)
library(ggprism)
library(patchwork)
library(meconetcomp)
library(magrittr)
library(igraph)
library(compositions)
library(microbiome)
library(phyloseq)
library(file2meco)
library(ggClusterNet)
library(WikidataR)
library(reshape2)
library(mecodev)
library(ggpubr)
library(factoextra)
library(FactoMineR)

```

# 2 Data import

## 2.1 Dataset generation

```{r}

b_otu <- read.csv("./b_OTU_V1.0.csv", header = TRUE, row.names = 1)
f_otu <- read.csv("./f_OTU_V1.0.csv", header = TRUE, row.names = 1)

sample <- read.csv("./envir_V1.0.csv", header = TRUE, row.names = 1)

b_tax <- read.csv("./b_tax_V1.0.csv", header = TRUE, row.names = 1)
f_tax <- read.csv("./f_tax_V1.0.csv", header = TRUE, row.names = 1)

```

```{r}

#Build a microtable and store it in a microtable object.
##bacteria
b_dataset <- microtable$new(sample_table = sample, 
                            otu_table = b_otu, 
                            tax_table = b_tax)
##fungi
f_dataset <- microtable$new(sample_table = sample, 
                            otu_table = f_otu, 
                            tax_table = f_tax)

b_dataset$tidy_dataset()
f_dataset$tidy_dataset()

b_dataset
f_dataset

```

## 2.2 Site dataset processing

```{r}

#Set the target site data to be analyzed
site_name <- "JFL"
#"GH"(boreal forest)、"DLS"(temperate forest)、"GNJ"(subtropical forest)、"JFL"(tropical forest)

```

### 2.2.1 Bacterial data of the target site

```{r}

b_dataset_filter <- clone(b_dataset)
b_dataset_filter$sample_table %<>% subset(Site == site_name)
b_dataset_filter$tidy_dataset()

#CK-bacteria
ck_b_dataset <- clone(b_dataset_filter)
ck_b_dataset$sample_table %<>% subset(Group == "CK")
ck_b_dataset$tidy_dataset()
ck_b_dataset$filter_taxa(freq = 0.75)

ck_b_dataset

#N50-bacteria
n50_b_dataset <- clone(b_dataset_filter)
n50_b_dataset$sample_table %<>% subset(Group == "N50")
n50_b_dataset$tidy_dataset()
n50_b_dataset$filter_taxa(freq = 0.75)

n50_b_dataset

#N100-bacteria
n100_b_dataset <- clone(b_dataset_filter)
n100_b_dataset$sample_table %<>% subset(Group == "N100")
n100_b_dataset$tidy_dataset()
n100_b_dataset$filter_taxa(freq = 0.75)
n100_b_dataset

```

### 2.2.2 Fungal data of the target site

```{r}

f_dataset_filter <- clone(f_dataset)
f_dataset_filter$sample_table %<>% subset(Site == site_name)
f_dataset_filter$tidy_dataset()

#CK-fungi
ck_f_dataset <- clone(f_dataset_filter)
ck_f_dataset$sample_table %<>% subset(Group == "CK")
ck_f_dataset$tidy_dataset()
ck_f_dataset$filter_taxa(freq = 0.5)
ck_f_dataset

#N50-fungi
n50_f_dataset <- clone(f_dataset_filter)
n50_f_dataset$sample_table %<>% subset(Group == "N50")
n50_f_dataset$tidy_dataset()
n50_f_dataset$filter_taxa(freq = 0.5)
n50_f_dataset

#N100-fungi
n100_f_dataset <- clone(f_dataset_filter)
n100_f_dataset$sample_table %<>% subset(Group == "N100")
n100_f_dataset$tidy_dataset()
n100_f_dataset$filter_taxa(freq = 0.5)
n100_f_dataset

```

```{r}

b_dataset_list <- list(CK = ck_b_dataset, N50 = n50_b_dataset, N100 = n100_b_dataset)
f_dataset_list <- list(CK = ck_f_dataset, N50 = n50_f_dataset, N100 = n100_f_dataset)

```

```{r}

#Set the taxa to be analyzed
group_name <- "f"# b = bacteria、f = fungi

if (group_name == "b") {
  data_list <- b_dataset_list
} else if (group_name == "f") {
  data_list <- f_dataset_list
} else if (group_name == "bf") {
  data_list <- bf_dataset_list
} else {
  stop("Invalid file_name value. Please use 'b', 'f', or 'bf'.")
}

file_name <- paste(site_name, group_name, sep = "_")

```

# 3 Construction of soil microbial network

```{r}

cor_ck <- trans_network$new(cor_method = NULL, dataset = data_list$CK)
cor_n50 <- trans_network$new(cor_method = NULL, dataset = data_list$N50)
cor_n100 <- trans_network$new(cor_method = NULL, dataset = data_list$N100)

```

```{r}

#CK
network_ck <- clone(cor_ck)
network_ck$cal_network(network_method = "SpiecEasi", 
                       SpiecEasi_method = "mb", 
                       lambda.min.ratio = 1e-2,
                       nlambda = 20
)

#N50
network_n50 <- clone(cor_n50)
network_n50$cal_network(network_method = "SpiecEasi", 
                        SpiecEasi_method = "mb", 
                        lambda.min.ratio = 1e-2,
                        nlambda = 20
)

#N100
network_n100 <- clone(cor_n100)
network_n100$cal_network(network_method = "SpiecEasi", 
                         SpiecEasi_method = "mb", 
                         lambda.min.ratio = 1e-2,
                         nlambda = 20
)

```


```{r}

#Network modular computing
network_ck$cal_module(method = "cluster_fast_greedy")
network_n50$cal_module(method = "cluster_fast_greedy")
network_n100$cal_module(method = "cluster_fast_greedy")

#Calculate adjacency matrix
network_ck$get_adjacency_matrix()
network_n50$get_adjacency_matrix()
network_n100$get_adjacency_matrix()

#Get section properties
network_ck$get_node_table(node_roles = TRUE)
network_n50$get_node_table(node_roles = TRUE)
network_n100$get_node_table(node_roles = TRUE)

#Get edge attributes
network_ck$get_edge_table()
network_n50$get_edge_table()
network_n100$get_edge_table()

```

```{r}

#Set file name
network_ck_file_name <- sprintf("%s_network_ck.gexf", file_name)
network_n50_file_name <- sprintf("%s_network_n50.gexf", file_name)
network_n100_file_name <- sprintf("%s_network_n100.gexf", file_name)

#Save network files and use them for visualization
network_ck$save_network(filepath = network_ck_file_name)
network_n50$save_network(filepath = network_n50_file_name)
network_n100$save_network(filepath = network_n100_file_name)

```

# 4 Calculation of network topology index

```{r}

network_list <- list()

network_list$CK <- network_ck
network_list$N50 <- network_n50
network_list$N100 <- network_n100

```

### 4.1 overall network properties

```{r}

graph_ck <- graph_from_data_frame(network_ck$res_edge_table, 
                                  vertices = network_ck$res_node_table, 
                                  directed = FALSE)

graph_n50 <- graph_from_data_frame(network_n50$res_edge_table, 
                                   vertices = network_n50$res_node_table, 
                                   directed = FALSE)

graph_n100 <- graph_from_data_frame(network_n100$res_edge_table, 
                                    vertices = network_n100$res_node_table, 
                                    directed = FALSE)

net_properties_ck <- net_properties.4(graph_ck)
net_properties_n50 <- net_properties.4(graph_n50)
net_properties_n100 <- net_properties.4(graph_n100)

net_properties_ck
net_properties_n50
net_properties_n100

```

```{r}

neg_edges_ck <- table(network_ck$res_edge_table$label)
neg_edges_n50 <- table(network_n50$res_edge_table$label)
neg_edges_n100 <- table(network_n100$res_edge_table$label)

calculate_neg_edge_ratio <- function(neg_edges_table) {
  num_neg_edges <- neg_edges_table["-"]
  num_pos_edges <- neg_edges_table["+"]
  total_edges <- sum(neg_edges_table)
  neg_edge_ratio <- num_neg_edges / total_edges
  return(neg_edge_ratio)
}

neg_edge_ratio_ck <- calculate_neg_edge_ratio(neg_edges_ck)
neg_edge_ratio_n50 <- calculate_neg_edge_ratio(neg_edges_n50)
neg_edge_ratio_n100 <- calculate_neg_edge_ratio(neg_edges_n100)

neg_edge_ratios <- data.frame(
  Group = c("CK", "N50", "N100"),
  Neg_Edge_Ratio = c(neg_edge_ratio_ck, neg_edge_ratio_n50, neg_edge_ratio_n100)
)

print(neg_edge_ratios)

```

### 4.2 Sample network properties

```{r}

np_sample_res <- subnet_property(network_list)
np_sample_res
colnames(np_sample_res)[colnames(np_sample_res) == "Sample"] <- "SampleID"
colnames(np_sample_res)[colnames(np_sample_res) == "Network"] <- "Group"

```


### 4.3 Cohesion calculation

```{r}

cohesion <- cohesionclass$new(network_list)
cohesion_microeco_result <- cohesion$res_list$sample
cohesion$plot(measure = "c_neg")
cohesion$plot(measure = "c_pos")

names(cohesion_microeco_result)[names(cohesion_microeco_result) == "sample"] <- "SampleID"

final_net_res <-  merge(np_sample_res, cohesion_microeco_result[, c("SampleID", "c_pos","c_neg")], by = "SampleID")
final_net_res$cohersion <- final_net_res$c_pos + abs(final_net_res$c_neg)

net_properties_file_name <- sprintf("%s_network_properties_SE.csv", file_name)
write.csv(final_net_res, file = net_properties_file_name)

colnames(final_net_res)

```






