R/bn_tabu_gen_1.R

In GmicR: Combines WGCNA and xCell readouts with bayesian network learrning to generate a Gene-Module Immune-Cell network (GMIC)

    #' Uses tabu search algorithm to learn the structure of discretized data.
    #' @import bnlearn
    #' @param Auto_WGCNA_OUTPUT an R object generated by Auto_WGCNA and 
    #' discretized using the Data_Prep function.
    #' @param score character string indicating the score used for structure 
    #' learning. If "bde" (default), prior is set to "uniform".  
    #' If bds is used, the prior is set to "marginal".
    #' @param bootstraps_replicates an integer for the number of 
    #' bootstraps_replicates
    #' used for structure learning. Default value is 500
    #' @param iss the imaginary sample size, used by the Bayesian 
    #' Dirichlet scores (bde and bds)
    #' It is also known as “equivalent sample size”.
    #' The default value is equal to 10.
    #' @examples 
    #' GMIC_Builder_disc_dir<-system.file("extdata", "GMIC_Builder_disc.Rdata", 
    #' package = "GmicR", mustWork = TRUE)
    #' load(GMIC_Builder_disc_dir)
    #' 
    #' no_cores<-1
    #' cl<-parallel::makeCluster(no_cores)
    #' 
    #' GMIC_net<-bn_tabu_gen(GMIC_Builder_disc,
    #' cluster = cl,
    #' bootstraps_replicates = 50, score = "bds")
    #' parallel::stopCluster(cl)
    #' @inheritParams bnlearn::boot.strength
    #' @inheritParams bnlearn::tabu
    #' @return The learned bayesian network
    #' @seealso \code{\link[bnlearn]{arc.strength}}
    #' @seealso \code{\link[bnlearn]{hc}}
    #' @seealso \code{\link[bnlearn]{score}}
    #' @export
    
    
    bn_tabu_gen<-function(Auto_WGCNA_OUTPUT,
    whitelist=NULL, blacklist=NULL,
    score = "bde", tabu = 50, iss = 10,
    cluster=NULL,
    debug = TRUE, bootstraps_replicates = 500){
    
    data<-Auto_WGCNA_OUTPUT$disc_data
    
    # mbde network start ------------------------------------------------------
    if(score == "bds"){
    prior = "marginal"} else if(score == "bde"){
    prior = "uniform"}
    
    tabu_net_intscore <-boot.strength(data, R = bootstraps_replicates, 
    algorithm = "tabu",cluster=cluster,
    debug = debug, algorithm.args = list(score = score,
    whitelist = whitelist,
    blacklist = blacklist,
    iss = iss, tabu = tabu, prior = prior))
    
    # Output preparation ------------------------------------------------------
    
    tabu_net_ave = averaged.network(tabu_net_intscore)
    Arcs_tabu_net_ave = arcs(tabu_net_ave)
    tabu_net_fitted<-c()
    
    if(directed(tabu_net_ave)){
    tabu_net_fitted<-bn.fit(tabu_net_ave, data, method = "bayes")
    } else{
    message("Graph partially directed, can't detect inverse relationships")
    }
    
    Output_tabu_net<-list(disc_data=data,
    score=score,
    prior=prior,
    tabu=tabu,
    iss=iss,
    bootstraps_replicates=bootstraps_replicates,
    whitelist=whitelist,
    blacklist=blacklist,
    custom_strength=tabu_net_intscore,
    ave.boot=tabu_net_ave,
    Arcs_tabu_net_ave=Arcs_tabu_net_ave,
    net_fitted=tabu_net_fitted,
    InverseR_Output=c())
    
    Auto_WGCNA_OUTPUT$Output_tabu_net<-Output_tabu_net
    return(Auto_WGCNA_OUTPUT)
    }
    
    #' Generates a subgraph from query nodes
    #' @param bn_output R object output from bn_tabu_gen function
    #' @param Node_ids vector containing the nodes for subgraph generation.
    #' @param relationship_type the relationship to be returned for the 
    #' specified query nodes.
    #' The options are "mb", "nbr", "parents", "children". Default 
    #' setting is "nbr".
    #' @return a subgraph containg the selected nodes and relationships.
    #' @export
    
    Batch_Net<-function(bn_output, Node_ids, relationship_type = "nbr"){
    
    ave.boot<-bn_output$ave.boot
    InverseR_Output<-bn_output$InverseR_Output
    
    # sets up selection criteria
    net_var<-data.frame(cbind(seq(1,4)), row.names = c("mb", "nbr", "parents", 
    "children"))
    net_var
    x<-net_var[relationship_type,]
    
    # makes list object for outout
    
    Net_lists <- setNames(vector(length(Node_ids), mode="list"), Node_ids)
    Net_lists
    
    # returns all nodes for networks
    for (i in Node_ids) {
    Net_lists[[i]]<-ave.boot$nodes[[i]][[x]]
    }
    Net_list_final<-unique(c(unlist(Net_lists), names(Net_lists)))
    sub_graph<-subgraph(ave.boot, Net_list_final) # makes subgraph
    
    # matching inverser arcs
    primary_arc_ids<-paste(InverseR_Output$from,InverseR_Output$to,sep = "_")
    primary_arc_ids
    
    sub_arcs_ids<-data.frame(arcs(sub_graph), stringsAsFactors = FALSE)
    sub_arcs_ids$code<-paste(sub_arcs_ids$from,sub_arcs_ids$to,sep = "_")
    
    sub_inverseR<-as.matrix(sub_arcs_ids[sub_arcs_ids$code %in% primary_arc_ids,
    seq(1,2)])
    
    
    Output_sub_graph<-list(ave.boot=sub_graph,
    InverseR_Output=sub_inverseR)
    
    
    return(Output_sub_graph)
    
    }
    
    #' Identifies arcs between nodes with inverse relationships
    #' @export
    #' @importFrom gRain querygrain
    #' @importFrom gRbase compile
    #' @param Output a data frame containing the output of 
    #' BN_Conditions function.
    #' @param threshold number indicating the maximum slope for 
    #' defining negative 
    #' relationships.
    #' Default level is -0.3.
    #' @examples GMIC_net_dir<-system.file("extdata", "GMIC_net.Rdata", 
    #' package = "GmicR", mustWork = TRUE)
    #' load(GMIC_net_dir)
    #' GMIC_Final<-InverseARCs(GMIC_net, threshold = -0.3)
    #' @return arcs with inverse relationships
    
    InverseARCs<-function(Output, threshold=-0.3){
    

    # checking version of bnlearn ---------------------------------------------
    bnlearn_version<-utils::packageVersion("bnlearn")
    
    if(bnlearn_version < "4.6"){
    message("Please update to bnlearn 4.6 from https://www.bnlearn.com/")
    }else if(bnlearn_version >= "4.6"){
    fitted<-Output$Output_tabu_net$net_fitted
    input_arcs<-arcs(fitted)
    
    # arcs to test ------------------------------------------------------------
    test_arcs<-data.frame(input_arcs, stringsAsFactors = FALSE)
    test_arcs$slope<-c(0)
    parent_list<-unique(test_arcs$from)
    for (i in parent_list) {
    
    # LOW
    q_L<-list("L")
    names(q_L)<-i
    j_L<-gRbase::compile(as.grain(mutilated(fitted, evidence = q_L)))
    
    # HIGH
    q_H<-list("H")
    names(q_H)<-i
    j_H<-gRbase::compile(as.grain(mutilated(fitted, evidence = q_H)))
    
    query_arcs<-test_arcs[test_arcs$from %in% i,]
    row_ids<-rownames(query_arcs)
    
    query_nodes<-c(i, query_arcs$to)
    
    LOW<-data.frame(gRain::querygrain(j_L, nodes = query_nodes[-1]))
    HIGH<-data.frame(gRain::querygrain(j_H, nodes = query_nodes[-1]))
    
    q_LOW<-data.frame(gRain::querygrain(j_L, nodes = query_nodes[1]))
    q_HIGH<-data.frame(gRain::querygrain(j_H, nodes = query_nodes[1]))
    
    all<-rbind(LOW, HIGH)
    q_ll<-rbind(q_LOW, q_HIGH)
    
    totaL<-na.omit(merge(all, q_ll, by= "row.names", all = TRUE))
    totaL$Row.names<-NULL
    
    cor_all<-cor(totaL)
    cor_all
    
    for(arc_id in row_ids){
    fr<-test_arcs[arc_id,]$from
    tot<-test_arcs[arc_id,]$to
    test_arcs[arc_id,]$slope<-cor_all[fr,tot]
    Inverse_relationships_arcs<-test_arcs[test_arcs$slope<threshold,seq(1,2)]
    }}
    
    print(test_arcs)
    Output$Output_tabu_net$InverseR_Output<-Inverse_relationships_arcs
    return(Output)
    }
    }