#' @title Identifying Differential Bacteria Based on a Network Workflow.
#' @description This is a function to identify differential bacteria between case and control data set.
#'
#' @param case_dir string.The directory of diseased data set or a single diseased file.
#' @param control_dir string.The directory of healthy data set or a single healthy file.
#' @param net_case_dir string.The directory of network correlation of diseased data set or a single diseased network file.
#' @param net_control_dir string.The directory of network correlation of healthy data set or a single healthy network file.
#' @param scaled logical.If TURE then the NetMoss score are scaled.
#' @param deepSplit numerical. For method "hybrid", can be either logical or integer in the range 0 to 4. For method "tree", must be logical. In both cases, provides a rough control over sensitivity to cluster splitting. The higher the value (or if TRUE), the more and smaller clusters will be produced. For the "hybrid" method, a finer control can be achieved via maxCoreScatter and minGap below.
#' @param minModuleSize numerical. Minimum module size. Set to 20 by default.
#' @param soft_pow numerical. Soft power for module division. Set to 5 by default.
#'
#' @return NetMoss score and module division results of each bacterium.
#' @export
#'
#'
NetMoss <-
  function(case_dir,
           control_dir,
           net_case_dir,
           net_control_dir,
           scaled = TRUE,
           deepSplit = 4,
           minModuleSize = 20,
           soft_pow = 5) {
		  
######reads files		  
getNetwork <-
  function(case_dir,
           control_dir,
           net_case_dir,
           net_control_dir) {
    print ("importing datasets...")

    ######net_case_dir is a data.frame#######
    if (length(net_case_dir) > 1 || length(net_control_dir) > 1)
    {
      ###read
      case_union = net_case_dir
      control_union = net_control_dir

      ##union
      union_genus1 <- rownames(case_union)
      union_genus2 <- rownames(control_union)

      union_genus <- union(union_genus1, union_genus2)

      union_matrix <-
        matrix(nrow = length(union_genus),
               ncol = length(union_genus),
               0)

      rownames(union_matrix) <- union_genus
      colnames(union_matrix) <- union_genus
      diag(union_matrix) <- 1

      ##combine
      #case
      fi <- case_union
      re_union <- union_matrix
      diff_genus <- setdiff(union_genus, rownames(fi))
      rownames(re_union) <- c(rownames(fi), diff_genus)
      colnames(re_union) <- c(rownames(fi), diff_genus)
      for (j in 1:nrow(fi)) {
        re_union[j, 1:nrow(fi)] <- as.numeric(fi[j,])
      }
      re_union <- re_union[union_genus, union_genus]
      diag(re_union) <- 1
      case_union_data_list <- re_union

      #control
      fi <- control_union
      re_union <- union_matrix
      diff_genus <- setdiff(union_genus, rownames(fi))
      rownames(re_union) <- c(rownames(fi), diff_genus)
      colnames(re_union) <- c(rownames(fi), diff_genus)
      for (j in 1:nrow(fi)) {
        re_union[j, 1:nrow(fi)] <- as.numeric(fi[j,])
      }
      re_union <- re_union[union_genus, union_genus]
      diag(re_union) <- 1
      control_union_data_list <- re_union

      integratedNet = list(case_union_data_list, control_union_data_list)

      return(integratedNet)
    } else if (file_test('-f', net_case_dir) &&
               file_test('-f', net_control_dir))
      #####net_case_dir is a file#####
    {
      ###read
      case_union = read.table(
        net_case_dir,
        header = T,
        sep = '\t',
        row.names = 1
      )
      control_union = read.table(
        net_control_dir,
        header = T,
        sep = '\t',
        row.names = 1
      )

      ##union
      union_genus1 <- rownames(case_union)
      union_genus2 <- rownames(control_union)

      union_genus <- union(union_genus1, union_genus2)

      union_matrix <-
        matrix(nrow = length(union_genus),
               ncol = length(union_genus),
               0)

      rownames(union_matrix) <- union_genus
      colnames(union_matrix) <- union_genus
      diag(union_matrix) <- 1

      ##combine
      #case
      fi <- case_union
      re_union <- union_matrix
      diff_genus <- setdiff(union_genus, rownames(fi))
      rownames(re_union) <- c(rownames(fi), diff_genus)
      colnames(re_union) <- c(rownames(fi), diff_genus)
      for (j in 1:nrow(fi)) {
        re_union[j, 1:nrow(fi)] <- as.numeric(fi[j,])
      }
      re_union <- re_union[union_genus, union_genus]
      diag(re_union) <- 1
      case_union_data_list <- re_union

      #control
      fi <- control_union
      re_union <- union_matrix
      diff_genus <- setdiff(union_genus, rownames(fi))
      rownames(re_union) <- c(rownames(fi), diff_genus)
      colnames(re_union) <- c(rownames(fi), diff_genus)
      for (j in 1:nrow(fi)) {
        re_union[j, 1:nrow(fi)] <- as.numeric(fi[j,])
      }
      re_union <- re_union[union_genus, union_genus]
      diag(re_union) <- 1
      control_union_data_list <- re_union

      integratedNet = list(case_union_data_list, control_union_data_list)

      return(integratedNet)

    } else
      ######else  combination ############
    {
      ######################################1. read data#################################################
      #case network
      setwd(case_dir)
      case_samples <- c()
      dir1 = list.files()
      n1 = length(dir1)
      for (i in 1:n1)
      {
        x.case = read.table(file = dir1[i],
                            header = T,
                            sep = "\t")
        case_samples <-
          c(case_samples, length(colnames(x.case)) - 1)
      }

      #control data
      setwd(control_dir)
      control_samples <- c()
      dir2 = list.files()
      n2 = length(dir2)
      for (i in 1:n2)
      {
        x.control = read.table(file = dir2[i],
                               header = T,
                               sep = "\t")
        control_samples <-
          c(control_samples, length(colnames(x.control)) - 1)
      }

      #case network
      case_data_list <- list()
      setwd(net_case_dir)
      dir3 = list.files()
      n3 = length(dir3)
      for (m in 1:n3)
      {
        x.case.net = read.table(
          file = dir3[m],
          header = T,
          sep = "\t",
          row.names = 1
        )
        case_data_list[[m]] = x.case.net
      }

      #control network
      control_data_list <- list()
      setwd(net_control_dir)
      dir4 = list.files()
      n4 = length(dir4)
      for (m in 1:n4)
      {
        x.control.net = read.table(
          file = dir4[m],
          header = T,
          sep = "\t",
          row.names = 1
        )
        control_data_list[[m]] = x.control.net
      }

      #union genus
      union_genus1 <- c()
      for (k in 1:length(case_data_list))
      {
        union_genus1 <- union(union_genus1, rownames(case_data_list[[k]]))
      }

      union_genus2 <- c()
      for (k in 1:length(control_data_list))
      {
        union_genus2 <-
          union(union_genus2, rownames(control_data_list[[k]]))
      }

      union_genus <- union(union_genus1, union_genus2)

      union_matrix <-
        matrix(nrow = length(union_genus),
               ncol = length(union_genus),
               0)

      rownames(union_matrix) <- union_genus
      colnames(union_matrix) <- union_genus
      diag(union_matrix) <- 1


      ######################################2. network construction#################################################
      print ("constructing networks...")
      case_union_data_list <- list()
      for (i in n1) {
        fi <- case_data_list[[i]]
        re_union <- union_matrix
        diff_genus <- setdiff(union_genus, rownames(fi))
        rownames(re_union) <- c(rownames(fi), diff_genus)
        colnames(re_union) <- c(rownames(fi), diff_genus)
        for (j in 1:nrow(fi)) {
          re_union[j, 1:nrow(fi)] <- as.numeric(fi[j,])
        }
        re_union <- re_union[union_genus, union_genus]
        diag(re_union) <- 1
        case_union_data_list[[i]] <- re_union
      }

      control_union_data_list <- list()
      for (i in n2) {
        fi <- control_data_list[[i]]
        re_union <- union_matrix
        diff_genus <- setdiff(union_genus, rownames(fi))
        rownames(re_union) <- c(rownames(fi), diff_genus)
        colnames(re_union) <- c(rownames(fi), diff_genus)
        for (j in 1:nrow(fi)) {
          re_union[j, 1:nrow(fi)] <- as.numeric(fi[j,])
        }
        re_union <- re_union[union_genus, union_genus]
        diag(re_union) <- 1
        control_union_data_list[[i]] <- re_union
      }


      ######################################3. combination#################################################
      print ("integrating networks...")
      pool_union <- union_matrix
      for (i in (1:(nrow(pool_union) - 1))) {
        for (j in ((i + 1):nrow(pool_union))) {
          son <- 0
          mom <- 0
          for (k in n1) {
            v <-
              (1 - case_union_data_list[[k]][i, j] ^ 2) / (case_samples[k] - 1)
            w <- 1 / v
            son <- son + w * case_union_data_list[[k]][i, j]
            mom <- mom + w
          }
          pool_union[i, j] <- son / mom
          pool_union[j, i] <- pool_union[i, j]
        }
      }
      case_union <- pool_union

      pool_union <- union_matrix
      for (i in (1:(nrow(pool_union) - 1))) {
        for (j in ((i + 1):nrow(pool_union))) {
          son <- 0
          mom <- 0
          for (k in n2) {
            v <-
              (1 - control_union_data_list[[k]][i, j] ^ 2) / (control_samples[k] - 1)
            w <- 1 / v
            son <- son + w * control_union_data_list[[k]][i, j]
            mom <- mom + w
          }
          pool_union[i, j] <- son / mom
          pool_union[j, i] <- pool_union[i, j]
        }
      }
      control_union <- pool_union

      integratedNet = list(case_union, control_union)

      return(integratedNet)
    }

  }

#####divide modules
divModule <-
  function(case_union,
           control_union,
           deepSplit = 4,
           minModuleSize = 20,
           soft_pow = 5) {
    print ("dividing modules...")

    adj_control <- as.matrix(0.5 + 0.5 * control_union)
    adj_case <- as.matrix(0.5 + 0.5 * case_union)

    pow_adjacency <- function(mat, pow) {
      mat_soft <- mat
      for (i in 1:nrow(mat)) {
        for (j in 1:ncol(mat)) {
          mat_soft[i, j] <- mat[i, j] ^ pow
        }
      }
      return(mat_soft)
    }

    ##health network modularization
    for (i in 7) {
      soft_pow_control <- soft_pow
      adj_soft_control <-
        pow_adjacency(adj_control, soft_pow_control)
      TOM_control <- TOMsimilarity(adj_soft_control)
      dissTOM_control <- 1 - TOM_control
      # Call the hierarchical clustering function
      geneTree_control = hclust(as.dist(dissTOM_control), method = "average")

      minModuleSize = minModuleSize

      # Module identification using dynamic tree cut:
      dynamicMods_control = cutreeDynamic(
        dendro = geneTree_control,
        distM = dissTOM_control,
        deepSplit = deepSplit,
        pamRespectsDendro = FALSE,
        minClusterSize = minModuleSize
      )

      print(table(dynamicMods_control))
      dynamicColors_control = labels2colors(dynamicMods_control)
    }
    ##disease network modularization
    for (i in 7) {
      soft_pow_case <- soft_pow
      adj_soft_case <- pow_adjacency(adj_case, soft_pow_case)
      TOM_case <- TOMsimilarity(adj_soft_case)
      dissTOM_case <- 1 - TOM_case
      # Call the hierarchical clustering function
      geneTree_case = hclust(as.dist(dissTOM_case), method = "average")

      # We like large modules, so we set the minimum module size relatively high:
      minModuleSize = minModuleSize

      # Module identification using dynamic tree cut:
      dynamicMods_case = cutreeDynamic(
        dendro = geneTree_case,
        distM = dissTOM_case,
        deepSplit = deepSplit,
        pamRespectsDendro = FALSE,
        minClusterSize = minModuleSize
      )

      print(table(dynamicMods_case))
      dynamicColors_case = labels2colors(dynamicMods_case)
    }

    modDivision = list(dissTOM_case,
                       dissTOM_control,
                       dynamicMods_case,
                       dynamicMods_control)

    return(modDivision)

  }

#####calculate NetMoss score
NetzGO <-
  function(control_mat,
           case_mat,
           control_dist,
           case_dist,
           control_mod,
           case_mod,
           scaled = TRUE) {
    if (scaled) {
      print("distance matrix scaled")
      control_dist_scaled <-
        (control_dist - min(control_dist)) / (max(control_dist) - min(control_dist))
      case_dist_scaled <-
        (case_dist - min(case_dist)) / (max(case_dist) - min(case_dist))
      control_dist <- control_dist_scaled
      case_dist <- case_dist_scaled
    }

    dynamicMods_inter <-
      paste("control_mod_", control_mod, "-", "case_mod_", case_mod, sep = "")
    dynamicMods_inter_names <- names(table(dynamicMods_inter))

    taxon_names <- rownames(control_mat)
    nodes <-
      data.frame(taxon_names, control_mod, case_mod, dynamicMods_inter)
    nodes$taxon_names <- as.character(nodes$taxon_names)
    nodes$control.degree <- rowSums(control_mat) - 1
    nodes$case.degree <- rowSums(case_mat) - 1

    dynamicMods_control_names <-
      as.numeric(names(table(control_mod)))
    dynamicMods_case_names <- as.numeric(names(table(case_mod)))

    nodes$control.mod.degree <- NA
    nodes$case.mod.degree <- NA
    for (i in dynamicMods_control_names) {
      row_num <- which(control_mod == i)
      col_num <- which(control_mod == i)
      nodes$control.mod.degree[row_num] <-
        rowSums(control_mat[row_num, col_num]) - 1 - rowSums(control_mat[row_num,-col_num])
    }

    for (i in dynamicMods_case_names) {
      row_num <- which(case_mod == i)
      col_num <- which(case_mod == i)
      nodes$case.mod.degree[row_num] <-
        rowSums(case_mat[row_num, col_num]) - 1 - rowSums(case_mat[row_num,-col_num])
    }

    nodes$diff.mod.degree <-
      (nodes$case.mod.degree - min(nodes$case.mod.degree)) / (max(nodes$case.mod.degree) -
                                                                min(nodes$case.mod.degree)) - (nodes$control.mod.degree - min(nodes$control.mod.degree)) /
      (max(nodes$control.mod.degree) - min(nodes$control.mod.degree))
    nodes$diff.mod.degree.abs <-
      abs((nodes$case.mod.degree - min(nodes$case.mod.degree)) / (max(nodes$case.mod.degree) -
                                                                    min(nodes$case.mod.degree)) - (nodes$control.mod.degree - min(nodes$control.mod.degree)) /
            (
              max(nodes$control.mod.degree) - min(nodes$control.mod.degree)
            )
      )

    control.nodes.to.mods.dist <-
      as.data.frame(matrix(
        nrow = length(taxon_names),
        ncol = length(dynamicMods_inter_names),
        NA
      ))
    rownames(control.nodes.to.mods.dist) <- taxon_names
    colnames(control.nodes.to.mods.dist) <- dynamicMods_inter_names
    case.nodes.to.mods.dist <- control.nodes.to.mods.dist

    nodes$netzgo <- 0
    for (i in dynamicMods_inter_names) {
      row_num <- which(dynamicMods_inter == i)
      for (j in row_num) {
        for (l in dynamicMods_inter_names) {
          col_num <- which(dynamicMods_inter == l)
          control.nodes.to.mods.dist[j, l] <-
            mean(control_dist[j, col_num])
          case.nodes.to.mods.dist[j, l] <-
            mean(case_dist[j, col_num])
        }
      }
    }
    diff.nodes.to.modes.dist <-
      case.nodes.to.mods.dist - control.nodes.to.mods.dist

    for (i in dynamicMods_inter_names) {
      target_num <- which(dynamicMods_inter == i)
      target_control_mod <- strsplit(i, "-")[[1]][1]
      target_case_mod <- strsplit(i, "-")[[1]][2]

      case_mod_names <-
        dynamicMods_inter_names[grep(target_control_mod, dynamicMods_inter_names)]
      case_mod_names <- setdiff(case_mod_names, i)
      control_mod_names <-
        dynamicMods_inter_names[grep(target_case_mod, dynamicMods_inter_names)]
      control_mod_names <- setdiff(control_mod_names, i)

      for (j in case_mod_names) {
        nodes$netzgo[target_num] <-
          nodes$netzgo[target_num] + diff.nodes.to.modes.dist[target_num, dynamicMods_inter_names %in% j]
      }
      for (k in control_mod_names) {
        nodes$netzgo[target_num] <-
          nodes$netzgo[target_num] - diff.nodes.to.modes.dist[target_num, dynamicMods_inter_names %in% k]
      }
    }
    nodes$netzgo.scale <- NA
    for (i in names(table(dynamicMods_inter))) {
      target_num <- which(nodes$dynamicMods_inter == i)
      for (j in target_num) {
        nodes$netzgo.scale[j] <-
          (nodes$netzgo[j] - min(nodes$netzgo[target_num])) / (max(nodes$netzgo[target_num]) - min(nodes$netzgo[target_num]))
      }
    }
    if (length(which(is.na(nodes$netzgo.scale))) != 0)
    {
      nodes = nodes[-which(is.na(nodes$netzgo.scale)), ]
    }


    nodes = nodes[, c("taxon_names",
                      "control_mod",
                      "case_mod",
                      "netzgo.scale")]
    colnames(nodes) = c("taxon_names",
                        "control_mod",
                        "case_mod",
                        "NetMoss_Score")

    #####p value
    p.table = data.frame()
    for (i in 1:nrow(case_mat))
    {
      ttax = rownames(case_mat)[i]
      p.table[i, 1] = ttax
      p.table[i, 2] = wilcox.test(case_dist[i, ], control_dist[i, ])$p.value
    }

    colnames(p.table) = c("taxon_names", "pvalue")
    rownames(nodes) = nodes$taxon_names
    p.table$score = nodes[as.character(p.table$taxon_names), "NetMoss_Score"]

    p.table$p.adj = p.adjust(p.table$pvalue, method = "BH")
    rownames(p.table) = p.table$taxon_names

    nodes$p.val = p.table[as.character(nodes$taxon_names), "pvalue"]
    nodes$p.adj = p.table[as.character(nodes$taxon_names), "p.adj"]

    nodes = nodes[which(nodes$p.adj < 0.05),]
    nodes = nodes[which(nodes$NetMoss_Score > 0),]

    return(nodes)

  }  


#####NetMoss		  
    my.wd = getwd()

    result = list()

    netAll = getNetwork(
      case_dir = case_dir,
      control_dir = control_dir,
      net_case_dir = net_case_dir,
      net_control_dir = net_control_dir
    )

    case_union = netAll[[1]]
    control_union = netAll[[2]]

    modAll = divModule(case_union = case_union,
                       control_union = control_union,
                       minModuleSize = minModuleSize,
                       deepSplit = deepSplit,
                       soft_pow = soft_pow)
    dissTOM_case = modAll[[1]]
    dissTOM_control = modAll[[2]]
    dynamicMods_case = modAll[[3]]
    dynamicMods_control = modAll[[4]]

    nodes_result = NetzGO(
      control_mat = control_union,
      case_mat = case_union,
      control_dist = dissTOM_control,
      case_dist = dissTOM_case,
      control_mod = as.numeric(dynamicMods_control),
      case_mod = as.numeric(dynamicMods_case),
      scaled = scaled
    )

    setwd(my.wd)

    result[[1]] = nodes_result
    result[[2]] = case_union
    result[[3]] = control_union

    return(result)

  }