R/net_node_cmp.R

In Guan06/mina: Microbial community dIversity and Network Analysis

###############################################################################

#' Compare the node features between networks.
#'
#' @param x The folder with all network inference results generated by bs_pm()
#' @param cmp The compared feature of node, default `contrast`.
#' @param dir The directory to store the alculated node features.
#' @examples
#' \dontrun{
#' net_node_cmp("./individual_bs_pm/", f = "contrast", dir = "./")
#'}
#' @rdname net_node_cmp
#' @export

setMethod("net_node_cmp", signature("character", "character"),
          function(x, cmp = "contrast", dir = "./"){
    
    if (!dir.exists(dir)) dir.create(dir)

    bs1_files <- sort(list.files(x, pattern = "_bs1.rds",
                                 full.names = TRUE))
    bs2_files <- sort(list.files(x, pattern = "_bs2.rds",
                                 full.names = TRUE))
    pm1_files <- sort(list.files(x, pattern = "_pm1.rds",
                                 full.names = TRUE))
    pm2_files <- sort(list.files(x, pattern = "_pm2.rds",
                                full.names = TRUE))
    len <- length(bs1_files)

    for (i in 1:len) {
        bs1 <- readRDS(bs1_files[i])
        bs2 <- readRDS(bs2_files[i])
        group_mn <- strsplit(basename(bs1_files[i]), "_bs1.rds")[[1]][1]
        group_m <- strsplit(group_mn, "_vs_")[[1]][1]
        group_n <- strsplit(group_mn, "_vs_")[[1]][2]

        y_bs <- list()
        y_bs[1] <- bs1
        y_bs[2] <- bs2

        this_m <- y_bs[[1]]
        this_n <- y_bs[[2]]

        ## calculate bootstrap distance
        bs_len <- length(this_m)
        node_feature_mn <- c()
        for (j1 in 1 : bs_len) {
            adj_m <- unlist(this_m[[j1]])
            adj_m[is.na(adj_m)] <- 0

            for (j2 in 1 : bs_len) {
                adj_n <- unlist(this_n[[j2]])
                adj_n[is.na(adj_n)] <- 0

                if (cmp == "contrast") {
                    this <- get_contrast(adj_m, adj_n)
                } else if (cmp == "ja") {
                    this <- get_ja(adj_m, adj_n)
                } else if (cmp == "ja0") {
                    this <- get_ja0(adj_m, adj_n)
                }

                node_feature_mn <- rbind(node_feature_mn, this)
            }
            colnames(node_feature_mn) <- rownames(adj_m)
        }
        dis_bs_file <- paste0(dir, "/dis_bs_", cmp, "_",
                              group_mn, ".rds")
        saveRDS(node_feature_mn, dis_bs_file, compress = "xz")

        ## the permutation results
        pm1 <- readRDS(pm1_files[i])
        pm2 <- readRDS(pm2_files[i])

        y_pm <- list()
        y_pm[1] <- pm1
        y_pm[2] <- pm2

        this_mp <- y_pm[[1]]
        this_np <- y_pm[[2]]

        pm_len <- length(this_mp)

        node_feature_mnp <- c()
        for (k1 in 1 : pm_len) {
            adj_mp <- unlist(this_mp[[k1]])
            adj_mp[is.na(adj_mp)] <- 0

            for (k2 in 1 : pm_len) {
                adj_np <- unlist(this_np[[k2]])
                adj_np[is.na(adj_np)] <- 0

                if (cmp == "contrast") {
                    this <- get_contrast(adj_mp, adj_np)
                } else if (cmp == "ja") {
                    this <- get_ja(adj_mp, adj_np)
                } else if (cmp == "ja0") {
                    this <- get_ja0(adj_mp, adj_np)
                }

                node_feature_mnp <- rbind(node_feature_mnp, this)
            }
            colnames(node_feature_mnp) <- rownames(adj_mp)
        }
        dis_pm_file <- paste0(dir, "/dis_pm_", cmp, "_",
                              group_mn, ".rds")
        saveRDS(node_feature_mnp, dis_pm_file, compress = "xz")
    }
})

###############################################################################

#' Get the contrast between two networks / adjacency matrices.
#' @param x The network/adjacency matrix of one of the compared condition.
#' @param y The other network/adjacency matrix.
#' @return The contrast between two matrices.
#' @keywords internal

get_contrast <- function(x, y) {
    rowSums(abs(x - y))
}

###############################################################################

#' Get the Jaccard distance between two networks / adjacency matrices.
#' @param x The network/adjacency matrix of one of the compared condition.
#' @param y The other network/adjacency matrix.
#' @return The Jaccard distance between two matrices.
#' @keywords internal

get_ja <- function(x, y) {
    t1 <- abs(x - y)
    t2 <- pmax(x, y)
    this <- rowSums(t1) / rowSums(t2)
}

###############################################################################

#' Get the Jaccard0 distance between two networks / adjacency matrices.
#' @param x The network/adjacency matrix of one of the compared condition.
#' @param y The other network/adjacency matrix.
#' @return The Jaccard0 distance between two matrices.
#' @keywords internal

get_ja0 <- function(x, y){
    t1 <- abs(x - y)
    t2 <- pmax(abs(x), abs(y))
    this <- rowSums(t1) / rowSums(t2)
}