Load the following libraries and install the package:

suppressPackageStartupMessages(library(doParallel))
suppressPackageStartupMessages(library(igraph))
## Warning: package 'igraph' was built under R version 4.2.3
suppressPackageStartupMessages(library(ggplot2))
## Warning: package 'ggplot2' was built under R version 4.2.3
suppressPackageStartupMessages(library(ggpubr))
suppressPackageStartupMessages(library(dagbagM))
suppressPackageStartupMessages(library(dagbag))
load("example.RData")
p<- dim(Y.n)[2] # no. of nodes
colnames(Y.n)<- paste0("p_",1:p)

temp.boot<- dagbagM::hc_boot_parallel(Y=Y.n, n.boot=10, nodeType=rep("c",p), whiteList=NULL, blackList=NULL, standardize=TRUE, tol = 1e-6, maxStep = 1000, restart=10, seed = 1,  nodeShuffle=TRUE, numThread = 2,verbose = FALSE)

res.agg<- dagbag::score_shd(temp.boot, alpha=1, threshold=0.5)

adj<- res.agg$adj.matrix
rownames(adj)<- colnames(adj)<- colnames(Y.n)

obj<- graph_from_adjacency_matrix(adj, mode = "directed", weighted = TRUE)
is_dag(obj)
## [1] TRUE
V(obj)$color<- c(rep("pink",(p-1)),"#99FFFF")

pdf("net.pdf",width = 30, height = 30, useDingbats=T)
plot.igraph(obj,vertex.size=10, edge.length=5, edge.color="black", vertex.label.cex=4,vertex.color=V(obj)$color,layout=layout_as_tree(obj, mode="all", circular = F))
dev.off()
## quartz_off_screen 
##                 2

Calculate the degree of the nodes and get the hub nodes

Degree of every node = number of the parent nodes + number of child nodes
child<- apply(adj,1,sum)
parent<- apply(adj,2,sum)
degree<- parent + child
degree.table<- cbind.data.frame("Nodes"=colnames(adj),"No. of child"=child, "No. of parent"=parent, "total"=degree)

degree.table$total = (degree.table$`No. of parent` +  degree.table$`No. of child`)
s<- sort(degree.table$total, ind=T, decreasing = T)
degree.table.sort<- degree.table[s$ix,]

# dat<- cbind.data.frame("Nodes" = rep(degree.table.sort$Nodes, 3), "edges" = c(degree.table.sort$total, degree.table.sort$parent, degree.table.sort$child), "type" = c(rep("Total",p), rep("Parent",p), rep("Child",p)))
# dat$Nodes<- factor(dat$Nodes, levels=unique(dat$Nodes))

# ggplot(dat, aes(y=edges, x=Nodes)) + 
#   geom_bar(position="stack", stat="identity") + theme_bw() +
#   theme(axis.text.x = element_text(colour="black",size=15,angle=45,hjust=1,vjust=1,face="plain"),axis.title.x = element_text(colour="black",size=25,angle=0,hjust=.5,vjust=.5,face="plain"), axis.text.y = element_text(colour="black",size=15,angle=0,hjust=1,vjust=0,face="plain"),  axis.title.y = element_text(colour="black",size=25,angle=90,hjust=.5,vjust=.5,face="plain"),legend.title = element_text(colour="black",size=15), legend.text=element_text(colour="black",size=15), legend.direction = "horizontal", legend.position = "top") + theme(plot.margin = unit(c(1,1,1,1), "cm"))

Extracting modules from the network

Modules are closely connected markers with similar functional roles/activities
load("adj_retro.RData")
obj<- graph_from_adjacency_matrix(adj, mode = "undirected", weighted = TRUE)

formodule<-cluster_edge_betweenness(obj, weights = NULL, directed = FALSE, edge.betweenness = TRUE, merges = TRUE, bridges = TRUE, modularity = TRUE, membership = TRUE)

module<-communities(formodule)
size<- unlist(lapply(module,length))

mod<- adj[which(rownames(adj) %in% module[[1]]), which(colnames(adj) %in% module[[1]])]
sum(mod)
## [1] 20
dim(mod)
## [1] 19 19
p=dim(mod)[1]
obj<- graph_from_adjacency_matrix(mod, mode = "directed", weighted = TRUE)
V(obj)$color<- rep("pink",p)

pdf("mod1.pdf",width = 30, height = 30, useDingbats=T)
plot.igraph(obj,vertex.size=10, edge.length=5, edge.color="black", vertex.label.cex=2,vertex.color=V(obj)$color,layout=layout_as_tree(obj, mode="all", circular = F))
dev.off()
## quartz_off_screen 
##                 2

Build DAG to infer causal association between markers in the module and clinical outcome

load("data_retro.RData")
data<-data[which(rownames(data) %in% rownames(mod)),]
load("sen_cptac.Rdata")
datay1<-cbind(t(data),"y"=sen.cptac)

load("prior.RData")
sub.adjn1<- sub.adjn[which(rownames(sub.adjn) %in% rownames(data)),which(rownames(sub.adjn) %in% rownames(data))]

sub.adjn1.ord<- sub.adjn1[order(match(rownames(sub.adjn1),rownames(data))),order(match(colnames(sub.adjn1),rownames(data)))]

sub.adjn2.ord<- cbind(sub.adjn1.ord,rep(0,dim(sub.adjn1)[2]))
sub.adjn3.ord<- rbind(sub.adjn2.ord,rep(0,dim(datay1)[2]))
colnames(sub.adjn3.ord)[dim(datay1)[2]]<-rownames(sub.adjn3.ord)[dim(datay1)[2]]<- "y"
sub.adjn4.ord<- ifelse(sub.adjn3.ord==1,T,F)

p<- dim(datay1)[2]
res.hc.m<- dagbagM::hc_boot_parallel(Y=datay1, n.boot=100, nodeType=c(rep("c",p-1), "b"),  whiteList=sub.adjn4.ord, blackList=NULL, standardize=TRUE, tol = 1e-6, maxStep = 1000, restart=10, seed = 1,  nodeShuffle=TRUE, numThread = 2,verbose = FALSE)

adj=score_shd(res.hc.m,threshold=0.4, whitelist = sub.adjn4.ord)$adj.matrix
rownames(adj)<- colnames(adj)<- colnames(datay1)

#obj<- graph_from_adjacency_matrix(adj, mode = "directed", weighted = TRUE)
#V(obj)$color<- c(rep("pink",(p-1)),"#99FFFF")

#pdf("adj_mod1.pdf",width = 30, height = 30, useDingbats=T)
#plot.igraph(obj,vertex.size=10, edge.length=5, edge.color="black", vertex.label.cex=3,vertex.color=V(obj)$color,layout=layout_as_tree(obj, mode="all", circular = F))
#dev.off()

#boxplot(datay1[,colnames(datay1)=="SAMM50"] ~ datay1[,colnames(datay1)=="y"])
#boxplot(datay1[,colnames(datay1)=="IMMT"] ~ datay1[,colnames(datay1)=="y"])