R/net_grp_cmp.R

In Guan06/mina: Microbial community dIversity and Network Analysis

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

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

setMethod("net_grp_cmp", signature("character", "ANY", "ANY", "ANY"),
    function(x, cmp = "contrast", dir = "./", grp){

    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)
    log <- c()

    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)
        grp_feature_mn <- c()
        for (j1 in 1 : bs_len) {
            adj_m <- unlist(this_m[[j1]])
            adj_m[is.na(adj_m)] <- 0

            log <- rbind(log,
                         paste0(group_mn, " bs_", j1, ": ", nrow(adj_m)))

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

                adj_m[is.na(adj_m)] <- 0
                adj_n[is.na(adj_n)] <- 0

                if (cmp == "contrast") {
                    this <- get_contrast_grp(adj_m, adj_n, grp)
                } else if (cmp == "ja") {
                    this <- get_ja_grp(adj_m, adj_n, grp)
                } else if (cmp == "ja0") {
                    this <- get_ja0_grp(adj_m, adj_n, grp)
                }

                grp_feature_mn <- rbind(grp_feature_mn, this)
            }
            colnames(grp_feature_mn) <- names(this)
        }
        dis_bs_file <- paste0(dir, "/dis_bs_", cmp, "_",
                              group_mn, ".rds")
        saveRDS(grp_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)

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

            log <- rbind(log, paste0(group_mn, " pm_", k1,
                              ": ", nrow(adj_mp)))

            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_grp(adj_mp, adj_np, grp)
                } else if (cmp == "ja") {
                    this <- get_ja_grp(adj_mp, adj_np, grp)
                } else if (cmp == "ja0") {
                    this <- get_ja0_grp(adj_mp, adj_np, grp)
                }

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

    }
)

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

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

get_contrast_grp <- function(x, y, grp) {
    grp$Group[is.na(grp$Group)] <- "Unassigned"
    group <- grp$Group

    x <- x[match(grp$ID, rownames(x)), ]
    y <- y[match(grp$ID, rownames(y)), ]
    z <- abs(x - y)

    contrast <- rowsum(z, group)
    contrast <- rowSums(contrast)
}

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

#' 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.
#' @param grp The table with group information.
#' @return The Jaccard distance between two matrices.
#' @keywords internal

get_ja_grp <- function(x, y, grp) {
    x <- x[match(grp$ID, rownames(x)), ]
    y <- y[match(grp$ID, rownames(y)), ]

    t1 <- abs(x - y)
    t2 <- pmax(x, y)
    
    grp$Group[is.na(grp$Group)] <- "Unassigned"
    group <- grp$Group

    t11 <- rowsum(t1, group)
    t22 <- rowsum(t2, group)

    this <- rowSums(t11) / rowSums(t22)
}
###############################################################################

#' 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.
#' @param grp The table with group information.
#' @return The Jaccard0 distance between two matrices.
#' @keywords internal

get_ja0_grp <- function(x, y, grp){
    x <- x[match(grp$ID, rownames(x)), ]
    y <- y[match(grp$ID, rownames(y)), ]

    t1 <- abs(x - y)
    t2 <- pmax(abs(x), abs(y))

    grp$Group[is.na(grp$Group)] <- "Unassigned"
    group <- grp$Group

    t11 <- rowsum(t1, group)
    t22 <- rowsum(t2, group)

    this <- rowSums(t11) / rowSums(t22)
}