diff --git a/publication_figures/src/makefig3old.R b/publication_figures/src/makefig3old.R
new file mode 100755
index 0000000..3a00f29
--- /dev/null
+++ b/publication_figures/src/makefig3old.R
@@ -0,0 +1,68 @@
+library(cowplot)
+library(ggplot2)
+theme_set(theme_cowplot(font_size=9))
+lambda <- seq(0.3,1.0, length.out = 50)
+cat.plot <- c("uncorrected", "RIN", "multi-covariate", "PC")
+# select category to plot
+select_category <- function(category_name, pr_table){
+  pr_table <- pr_table[which(pr_table$type %in% category_name),]
+  pr_table$type <- factor(pr_table$type, levels = category_name)
+  pr_table$lambda <- rep(lambda,length(category_name))
+  pr_table
+}
+plot.thyroid <- readRDS("/work-zfs/abattle4/parsana/networks_correction/results/PR/pr_density_glasso_canonical_thyroid.Rds")
+plot.thyroid <- select_category(cat.plot, plot.thyroid)
+plot.thyroid <- ggplot(plot.thyroid, aes(x = lambda, y = 1-precision, colour = type)) + geom_point(size =0.3) +
+  xlab("lambda") + ylab("FDR")+ggtitle("Thyroid") 
+plot.lung <- readRDS("/work-zfs/abattle4/parsana/networks_correction/results/PR/pr_density_glasso_canonical_lung.Rds")
+plot.lung <- select_category(cat.plot, plot.lung)
+plot.lung <- ggplot(plot.lung, aes(x = lambda, y = 1-precision, colour = type)) + geom_point(size = 0.3) +
+  xlab("lambda") + ylab("FDR")+ggtitle("Lung")
+plot.blood <- readRDS("/work-zfs/abattle4/parsana/networks_correction/results/PR/pr_density_glasso_canonical_blood.Rds")
+plot.blood <- select_category(cat.plot, plot.blood)
+plot.blood <- ggplot(plot.blood, aes(x = lambda, y = 1-precision, colour = type)) + geom_point(size = 0.3) +
+  xlab("lambda") + ylab("FDR")+ggtitle("Whole Blood")
+plot.skin <- readRDS("/work-zfs/abattle4/parsana/networks_correction/results/PR/pr_density_glasso_canonical_skin.Rds")
+plot.skin <- select_category(cat.plot, plot.skin)
+plot.skin <- ggplot(plot.skin, aes(x = lambda, y = 1-precision, colour = type)) + geom_point(size = 0.3) +
+  xlab("lambda") + ylab("FDR")+ggtitle("skin")
+plot.muscle <- readRDS("/work-zfs/abattle4/parsana/networks_correction/results/PR/pr_density_glasso_canonical_muscle.Rds")
+plot.muscle <- select_category(cat.plot, plot.muscle)
+plot.muscle <- ggplot(plot.muscle, aes(x = lambda, y = 1-precision, colour = type)) + geom_point(size = 0.3) +
+  xlab("lambda") + ylab("FDR")+ggtitle("muscle")
+plot.subcutaneous <- readRDS("/work-zfs/abattle4/parsana/networks_correction/results/PR/pr_density_glasso_canonical_subcutaneous.Rds")
+plot.subcutaneous <- select_category(cat.plot, plot.subcutaneous)
+plot.subcutaneous <- ggplot(plot.subcutaneous, aes(x = lambda, y = 1-precision, colour = type)) + geom_point(size = 0.3) +
+  xlab("lambda") + ylab("FDR")+ggtitle("subcutaneous")
+plot.artery <- readRDS("/work-zfs/abattle4/parsana/networks_correction/results/PR/pr_density_glasso_canonical_artery.Rds")
+plot.artery <- select_category(cat.plot, plot.artery)
+plot.artery <- ggplot(plot.artery, aes(x = lambda, y = 1-precision, colour = type)) + geom_point(size = 0.3) +
+  xlab("lambda") + ylab("FDR")+ggtitle("artery")
+plot.nerve <- readRDS("/work-zfs/abattle4/parsana/networks_correction/results/PR/pr_density_glasso_canonical_nerve.Rds")
+plot.nerve <- select_category(cat.plot, plot.nerve)
+plot.nerve <- ggplot(plot.nerve, aes(x = lambda, y = 1-precision, colour = type)) + geom_point(size = 0.3) +
+  xlab("lambda") + ylab("FDR")+ggtitle("nerve")
+###############################
+fig3 <- plot_grid(plot.blood + xlim(0.3,1.0) + ylim(0, 1) + theme(legend.position="none"),
+                  plot.lung + xlim(0.3,1.0) + ylim(0, 1) + theme(legend.position="none"),
+                  # plot.muscle + xlim(0.3,1.0) + ylim(0, 1) + theme(legend.position="none"),
+                  # plot.artery + xlim(0.3,1.0) + ylim(0, 1) + theme(legend.position="none"),
+                  # plot.skin + xlim(0.3,1.0) + ylim(0, 1) + theme(legend.position="none"),
+                  # plot.nerve + xlim(0.3,1.0) + ylim(0, 1) + theme(legend.position="none"),
+                  # plot.subcutaneous + xlim(0.3,1.0) + ylim(0, 1) + theme(legend.position="none"),
+                  plot.thyroid + xlim(0.3,1.0) + ylim(0, 1) +  theme(legend.position="none"),
+                  align = 'vh',
+                  labels = c("a", "b", "c"),
+                  hjust = -1,
+                  nrow = 1
+)
+legend <- get_legend(plot.lung +
+                       theme(legend.key = element_rect(color = "black", linetype = "solid", size = 0.5),
+                             legend.key.size = unit(0.3, "cm"), legend.key.height=unit(1.5,"line")) +
+                       guides(colour = guide_legend(override.aes = list(size= 1))))
+fig3a <- plot_grid( fig3, legend, rel_widths = c(3, .4))
+
+pdf("/work-zfs/abattle4/parsana/networks_correction/publication_figures/fig3.pdf", height = 2.5, width = 7.2)
+print(fig3a)
+dev.off()
+
diff --git a/publication_figures/src/supp_fig13.R b/publication_figures/src/supp_fig13.R
new file mode 100644
index 0000000..180c686
--- /dev/null
+++ b/publication_figures/src/supp_fig13.R
@@ -0,0 +1,86 @@
+library(cowplot)
+library(ggplot2)
+library(parallel)
+library(igraph)
+library(reshape2)
+theme_set(theme_cowplot(font_size=9))
+
+create.igraph <- function(network){
+    # network[which(network != 0)] <- 1
+    diag(network) <- 0
+    network[network!=0] <- 1
+    graph_from_adjacency_matrix(network, mode = "undirected")
+}
+
+
+fn <- dir("../../networks/", recursive=T, full.names = T)
+fn <- fn[grep("glasso", fn)]
+fn <- fn[grep(paste("raw", "rin", "mc", "/pc/", sep = "|"), fn)]
+
+igraph_net <- mclapply(fn, function(x){
+        print(x)
+        load(x)
+        net.igraph <- lapply(dat.net[1:10], function(eachNet){
+                eachNet[lower.tri(eachNet)] <- t(eachNet)[lower.tri(eachNet)]
+                create.igraph(eachNet)
+                })
+        net.igraph
+        }, mc.cores = 16)
+
+# igraph_net <- lapply(igraph_net, function(x,y){
+#         y <- y
+#         this_net <- lapply(x, function(thisx, thisy){
+#                 V(thisx)$name <- thisy
+#                 thisx
+#                 }, y)
+#         this_net
+#         }, dat.gene.symbol$gene_symbol)
+
+
+nodes_with_20 <- mclapply(igraph_net, function(x)
+                sapply(x, function(y) sum(degree(y) >=20)), mc.cores = 16)
+
+names(nodes_with_20) <- fn
+
+raw <- nodes_with_20[grep("raw", fn)]
+rin <- nodes_with_20[grep("rin", fn)]
+mc <- nodes_with_20[grep("mc", fn)]
+pc <- nodes_with_20[grep("/pc/", fn)]
+
+hub_net <- melt(data.frame(raw, rin, mc, pc))
+
+hub_net$tissue <- NA
+hub_net$tissue[grep("Artery", hub_net$variable)] <- "Tibial Artery"
+hub_net$tissue[grep("Blood", hub_net$variable)] <- "Whole Blood"
+hub_net$tissue[grep("Lung", hub_net$variable)] <- "Lung"
+hub_net$tissue[grep("Muscle", hub_net$variable)] <- "Muscle"
+hub_net$tissue[grep("Subcutaneous", hub_net$variable)] <- "Adipose Subcutaneous"
+hub_net$tissue[grep("Nerve", hub_net$variable)] <- "Nerve Tibial"
+hub_net$tissue[grep("Thyroid", hub_net$variable)] <- "Thyroid"
+hub_net$tissue[grep("Skin", hub_net$variable)] <- "Sun exposed Skin"
+
+hub_net$version <- NA
+hub_net$version[grep("raw", hub_net$variable)] <- "Uncorrected"
+hub_net$version[grep("rin", hub_net$variable)] <- "RIN"
+hub_net$version[grep("mc", hub_net$variable)] <- "multi-covariate"
+hub_net$version[grep("..pc.", hub_net$variable)] <- "PC"
+
+
+select_category <- function(category_name, pr_table){
+  pr_table <- pr_table[which(pr_table$version %in% category_name),]
+  pr_table$version <- factor(pr_table$version, levels = category_name)
+  pr_table
+}
+
+which_version <- c("Uncorrected", "RIN", "multi-covariate", "PC")
+
+plot_this <- select_category(which_version, hub_net)
+
+thisplot <- ggplot(aes(y = value, x = tissue, fill = version),
+        data = plot_this)+
+geom_boxplot()+ ylab("# Nodes with >=20 neighbors")+ xlab("Tissue") + theme_classic()+theme(axis.text.x = element_text(angle = 90, hjust = 1))
+
+pdf("../suppfig13.pdf", height = 7, width = 12)
+print(thisplot)
+dev.off()
+
diff --git a/publication_figures/src/supp_fig14.R b/publication_figures/src/supp_fig14.R
new file mode 100644
index 0000000..a62b2fd
--- /dev/null
+++ b/publication_figures/src/supp_fig14.R
@@ -0,0 +1,136 @@
+library(cowplot)
+library(ggplot2)
+theme_set(theme_cowplot(font_size=9))
+cutheights <- seq(0.9,1.0, length.out = 50)
+
+cat.plot <- c("uncorrected", "RIN", "multi-covariate", "PC")
+
+# select category to plot
+select_category <- function(category_name, pr_table){
+  pr_table <- pr_table[which(pr_table$type %in% category_name),]
+  pr_table$type <- factor(pr_table$type, levels = category_name)
+  pr_table$cutheights <- rep(cutheights,length(category_name))
+  pr_table$fnr <- pr_table$false_negatives/ ((4978 * (4978 - 1))/2 - pr_table$density)
+  pr_table
+}
+
+plot.thyroid <- readRDS("/work-zfs/abattle4/parsana/networks_correction/results/PR/pr_density_wgcna-signed_canonical_thyroid.Rds")
+plot.thyroid <- select_category(cat.plot, plot.thyroid)
+plot.thyroid <- ggplot(plot.thyroid, aes(x = cutheights, y  = fnr, colour = type)) + geom_point(size =0.3) +
+  xlab("cut-heights") + ylab("FNR")+ggtitle("Thyroid")
+plot.lung <- readRDS("/work-zfs/abattle4/parsana/networks_correction/results/PR/pr_density_wgcna-signed_canonical_lung.Rds")
+plot.lung <- select_category(cat.plot, plot.lung)
+plot.lung <- ggplot(plot.lung, aes(x = cutheights, y  = fnr, colour = type)) + geom_point(size = 0.3) +
+  xlab("cut-heights") + ylab("FNR")+ggtitle("Lung")
+plot.blood <- readRDS("/work-zfs/abattle4/parsana/networks_correction/results/PR/pr_density_wgcna-signed_canonical_blood.Rds")
+plot.blood <- select_category(cat.plot, plot.blood)
+plot.blood <- ggplot(plot.blood, aes(x = cutheights, y  = fnr, colour = type)) + geom_point(size = 0.3) +
+  xlab("cut-heights") + ylab("FNR")+ggtitle("Whole Blood")
+plot.skin <- readRDS("/work-zfs/abattle4/parsana/networks_correction/results/PR/pr_density_wgcna-signed_canonical_skin.Rds")
+plot.skin <- select_category(cat.plot, plot.skin)
+plot.skin <- ggplot(plot.skin, aes(x = cutheights, y  = fnr, colour = type)) + geom_point(size = 0.3) +
+  xlab("cut-heights") + ylab("FNR")+ggtitle("skin")
+plot.muscle <- readRDS("/work-zfs/abattle4/parsana/networks_correction/results/PR/pr_density_wgcna-signed_canonical_muscle.Rds")
+plot.muscle <- select_category(cat.plot, plot.muscle)
+plot.muscle <- ggplot(plot.muscle, aes(x = cutheights, y  = fnr, colour = type)) + geom_point(size = 0.3) +
+  xlab("cut-heights") + ylab("FNR")+ggtitle("muscle")
+plot.subcutaneous <- readRDS("/work-zfs/abattle4/parsana/networks_correction/results/PR/pr_density_wgcna-signed_canonical_subcutaneous.Rds")
+plot.subcutaneous <- select_category(cat.plot, plot.subcutaneous)
+plot.subcutaneous <- ggplot(plot.subcutaneous, aes(x = cutheights, y  = fnr, colour = type)) + geom_point(size = 0.3) +
+  xlab("cut-heights") + ylab("FNR")+ggtitle("subcutaneous")
+plot.artery <- readRDS("/work-zfs/abattle4/parsana/networks_correction/results/PR/pr_density_wgcna-signed_canonical_artery.Rds")
+plot.artery <- select_category(cat.plot, plot.artery)
+plot.artery <- ggplot(plot.artery, aes(x = cutheights, y  = fnr, colour = type)) + geom_point(size = 0.3) +
+  xlab("cut-heights") + ylab("FNR")+ggtitle("artery")
+plot.nerve <- readRDS("/work-zfs/abattle4/parsana/networks_correction/results/PR/pr_density_wgcna-signed_canonical_nerve.Rds")
+plot.nerve <- select_category(cat.plot, plot.nerve)
+plot.nerve <- ggplot(plot.nerve, aes(x = cutheights, y  = fnr, colour = type)) + geom_point(size = 0.3) +
+  xlab("cut-heights") + ylab("FNR")+ggtitle("nerve")
+###############################
+fig2 <- plot_grid(plot.blood + xlim(0.9,1.0) + ylim(0.045, 0.050) + theme(legend.position="none"),
+                  plot.lung + xlim(0.9,1.0) + ylim(0.045, 0.050) + theme(legend.position="none"),
+                  # plot.muscle + xlim(0.9,1.0) + + theme(legend.position="none"),
+                  # plot.artery + xlim(0.9,1.0) + + theme(legend.position="none"),
+                  # plot.skin + xlim(0.9,1.0) + + theme(legend.position="none"),
+                  # plot.nerve + xlim(0.9,1.0) + + theme(legend.position="none"),
+                  # plot.subcutaneous + xlim(0.9,1.0) + + theme(legend.position="none"),
+                  plot.thyroid + xlim(0.9,1.0) + ylim(0.045, 0.050) + theme(legend.position="none"),
+                  align = 'vh',
+                  labels = c("a", "b", "c"),
+                  hjust = -1,
+                  nrow = 1
+)
+legend <- get_legend(plot.lung +
+                       theme(legend.key = element_rect(color = "black", linetype = "solid", size = 0.5),
+                             legend.key.size = unit(0.3, "cm"), legend.key.height=unit(1.5,"line")) +
+                       guides(colour = guide_legend(override.aes = list(size= 1))))
+fig2a <- plot_grid( fig2, legend, rel_widths = c(3, .4))
+
+
+####################################### Fig 2b glasso ###################################
+lambda <- seq(0.3,1.0, length.out = 50)
+cat.plot <- c("uncorrected", "RIN", "multi-covariate", "PC")
+# select category to plot
+select_category <- function(category_name, pr_table){
+  pr_table <- pr_table[which(pr_table$type %in% category_name),]
+  pr_table$type <- factor(pr_table$type, levels = category_name)
+  pr_table$lambda <- rep(lambda,length(category_name))
+  pr_table$fnr <- pr_table$false_negatives/ ((4978 * (4978 - 1))/2 - pr_table$density)
+  pr_table
+}
+plot.thyroid <- readRDS("/work-zfs/abattle4/parsana/networks_correction/results/PR/pr_density_glasso_canonical_thyroid.Rds")
+plot.thyroid <- select_category(cat.plot, plot.thyroid)
+plot.thyroid <- ggplot(plot.thyroid, aes(x = lambda, y  = fnr, colour = type)) + geom_point(size =0.3) +
+  xlab("lambda") + ylab("FNR")+ggtitle("Thyroid")
+plot.lung <- readRDS("/work-zfs/abattle4/parsana/networks_correction/results/PR/pr_density_glasso_canonical_lung.Rds")
+plot.lung <- select_category(cat.plot, plot.lung)
+plot.lung <- ggplot(plot.lung, aes(x = lambda, y  = fnr, colour = type)) + geom_point(size = 0.3) +
+  xlab("lambda") + ylab("FNR")+ggtitle("Lung")
+plot.blood <- readRDS("/work-zfs/abattle4/parsana/networks_correction/results/PR/pr_density_glasso_canonical_blood.Rds")
+plot.blood <- select_category(cat.plot, plot.blood)
+plot.blood <- ggplot(plot.blood, aes(x = lambda, y  = fnr, colour = type)) + geom_point(size = 0.3) +
+  xlab("lambda") + ylab("FNR")+ggtitle("Whole Blood")
+plot.skin <- readRDS("/work-zfs/abattle4/parsana/networks_correction/results/PR/pr_density_glasso_canonical_skin.Rds")
+plot.skin <- select_category(cat.plot, plot.skin)
+plot.skin <- ggplot(plot.skin, aes(x = lambda, y  = fnr, colour = type)) + geom_point(size = 0.3) +
+  xlab("lambda") + ylab("FNR")+ggtitle("skin")
+plot.muscle <- readRDS("/work-zfs/abattle4/parsana/networks_correction/results/PR/pr_density_glasso_canonical_muscle.Rds")
+plot.muscle <- select_category(cat.plot, plot.muscle)
+plot.muscle <- ggplot(plot.muscle, aes(x = lambda, y  = fnr, colour = type)) + geom_point(size = 0.3) +
+  xlab("lambda") + ylab("FNR")+ggtitle("muscle")
+plot.subcutaneous <- readRDS("/work-zfs/abattle4/parsana/networks_correction/results/PR/pr_density_glasso_canonical_subcutaneous.Rds")
+plot.subcutaneous <- select_category(cat.plot, plot.subcutaneous)
+plot.subcutaneous <- ggplot(plot.subcutaneous, aes(x = lambda, y  = fnr, colour = type)) + geom_point(size = 0.3) +
+  xlab("lambda") + ylab("FNR")+ggtitle("subcutaneous")
+plot.artery <- readRDS("/work-zfs/abattle4/parsana/networks_correction/results/PR/pr_density_glasso_canonical_artery.Rds")
+plot.artery <- select_category(cat.plot, plot.artery)
+plot.artery <- ggplot(plot.artery, aes(x = lambda, y  = fnr, colour = type)) + geom_point(size = 0.3) +
+  xlab("lambda") + ylab("FNR")+ggtitle("artery")
+plot.nerve <- readRDS("/work-zfs/abattle4/parsana/networks_correction/results/PR/pr_density_glasso_canonical_nerve.Rds")
+plot.nerve <- select_category(cat.plot, plot.nerve)
+plot.nerve <- ggplot(plot.nerve, aes(x = lambda, y  = fnr, colour = type)) + geom_point(size = 0.3) +
+  xlab("lambda") + ylab("FNR")+ggtitle("nerve")
+###############################
+fig3 <- plot_grid(plot.blood + xlim(0.3,1.0) + ylim(0.045, 0.050)+ theme(legend.position="none"),
+                  plot.lung + xlim(0.3,1.0) + ylim(0.045, 0.050)+ theme(legend.position="none"),
+                  # plot.muscle + xlim(0.3,1.0) + theme(legend.position="none"),
+                  # plot.artery + xlim(0.3,1.0) + theme(legend.position="none"),
+                  # plot.skin + xlim(0.3,1.0) + theme(legend.position="none"),
+                  # plot.nerve + xlim(0.3,1.0) + theme(legend.position="none"),
+                  # plot.subcutaneous + xlim(0.3,1.0) + theme(legend.position="none"),
+                  plot.thyroid + xlim(0.3,1.0) + ylim(0.045, 0.050) + theme(legend.position="none"),
+                  align = 'vh',
+                  labels = c("a", "b", "c"),
+                  hjust = -1,
+                  nrow = 1
+)
+
+legend <- get_legend(plot.lung +
+                       theme(legend.key = element_rect(color = "black", linetype = "solid", size = 0.5),
+                             legend.key.size = unit(0.3, "cm"), legend.key.height=unit(1.5,"line")) +
+                       guides(colour = guide_legend(override.aes = list(size= 1))))
+fig2b <- plot_grid( fig3, legend, rel_widths = c(3, .4))
+
+pdf("../suppfig14.pdf", height = 5, width = 7.2)
+plot_grid(fig2a, fig2b, nrow = 2)
+dev.off()
diff --git a/publication_figures/src/supp_fig8.R b/publication_figures/src/supp_fig8.R
new file mode 100755
index 0000000..45e1cb9
--- /dev/null
+++ b/publication_figures/src/supp_fig8.R
@@ -0,0 +1,71 @@
+library(cowplot)
+library(ggplot2)
+theme_set(theme_cowplot(font_size=9))
+cutheights <- seq(0.9,1.0, length.out = 50)
+cat.plot <- c("uncorrected", "RIN", "multi-covariate", "PC")
+
+# select category to plot
+select_category <- function(category_name, pr_table){
+  pr_table <- pr_table[which(pr_table$type %in% category_name),]
+  pr_table$type <- factor(pr_table$type, levels = category_name)
+  pr_table$cutheights <- rep(cutheights,length(category_name))
+  pr_table
+}
+plot.thyroid <- readRDS("/work-zfs/abattle4/parsana/networks_correction/results/PR/pr_density_wgcna-signed_shared_thyroid.Rds")
+plot.thyroid <- select_category(cat.plot, plot.thyroid)
+plot.thyroid <- ggplot(plot.thyroid, aes(x = cutheights, y = 1-precision, colour = type)) + geom_point(size =0.3) +
+  xlab("cut-heights") + ylab("FDR")+ggtitle("Thyroid") 
+plot.lung <- readRDS("/work-zfs/abattle4/parsana/networks_correction/results/PR/pr_density_wgcna-signed_shared_lung.Rds")
+plot.lung <- select_category(cat.plot, plot.lung)
+plot.lung <- ggplot(plot.lung, aes(x = cutheights, y = 1-precision, colour = type)) + geom_point(size = 0.3) +
+  xlab("cut-heights") + ylab("FDR")+ggtitle("Lung")
+plot.blood <- readRDS("/work-zfs/abattle4/parsana/networks_correction/results/PR/pr_density_wgcna-signed_shared_blood.Rds")
+plot.blood <- select_category(cat.plot, plot.blood)
+plot.blood <- ggplot(plot.blood, aes(x = cutheights, y = 1-precision, colour = type)) + geom_point(size = 0.3) +
+  xlab("cut-heights") + ylab("FDR")+ggtitle("Whole Blood")
+plot.skin <- readRDS("/work-zfs/abattle4/parsana/networks_correction/results/PR/pr_density_wgcna-signed_shared_skin.Rds")
+plot.skin <- select_category(cat.plot, plot.skin)
+plot.skin <- ggplot(plot.skin, aes(x = cutheights, y = 1-precision, colour = type)) + geom_point(size = 0.3) +
+  xlab("cut-heights") + ylab("FDR")+ggtitle("skin")
+plot.muscle <- readRDS("/work-zfs/abattle4/parsana/networks_correction/results/PR/pr_density_wgcna-signed_shared_muscle.Rds")
+plot.muscle <- select_category(cat.plot, plot.muscle)
+plot.muscle <- ggplot(plot.muscle, aes(x = cutheights, y = 1-precision, colour = type)) + geom_point(size = 0.3) +
+  xlab("cut-heights") + ylab("FDR")+ggtitle("muscle")
+plot.subcutaneous <- readRDS("/work-zfs/abattle4/parsana/networks_correction/results/PR/pr_density_wgcna-signed_shared_subcutaneous.Rds")
+plot.subcutaneous <- select_category(cat.plot, plot.subcutaneous)
+plot.subcutaneous <- ggplot(plot.subcutaneous, aes(x = cutheights, y = 1-precision, colour = type)) + geom_point(size = 0.3) +
+  xlab("cut-heights") + ylab("FDR")+ggtitle("subcutaneous")
+plot.artery <- readRDS("/work-zfs/abattle4/parsana/networks_correction/results/PR/pr_density_wgcna-signed_shared_artery.Rds")
+plot.artery <- select_category(cat.plot, plot.artery)
+plot.artery <- ggplot(plot.artery, aes(x = cutheights, y = 1-precision, colour = type)) + geom_point(size = 0.3) +
+  xlab("cut-heights") + ylab("FDR")+ggtitle("artery")
+plot.nerve <- readRDS("/work-zfs/abattle4/parsana/networks_correction/results/PR/pr_density_wgcna-signed_shared_nerve.Rds")
+plot.nerve <- select_category(cat.plot, plot.nerve)
+plot.nerve <- ggplot(plot.nerve, aes(x = cutheights, y = 1-precision, colour = type)) + geom_point(size = 0.3) +
+  xlab("cut-heights") + ylab("FDR")+ggtitle("nerve")
+###############################
+legend <- get_legend(plot.lung +
+                       theme(legend.key = element_rect(color = "black", linetype = "solid", size = 0.5),
+                             legend.key.size = unit(0.3, "cm"), legend.key.height=unit(1.5,"line")) +
+                       guides(colour = guide_legend(override.aes = list(size= 1))))
+
+
+fig2 <- plot_grid(plot.blood + xlim(0.9,1.0) + ylim(0.5, 1) + theme(legend.position="none"),
+                  plot.lung + xlim(0.9,1.0) + ylim(0.5, 1) + theme(legend.position="none"),
+                  plot.thyroid + xlim(0.9,1.0) + ylim(0.5, 1) + theme(legend.position="none"),
+                  plot.muscle + xlim(0.9,1.0) + ylim(0.5, 1) + theme(legend.position="none"),
+                  plot.artery + xlim(0.9,1.0) + ylim(0.5, 1) + theme(legend.position="none"),
+                  plot.skin + xlim(0.9,1.0) + ylim(0.5, 1) + theme(legend.position="none"),
+                  plot.nerve + xlim(0.9,1.0) + ylim(0.5, 1) + theme(legend.position="none"),
+                  plot.subcutaneous + xlim(0.9,1.0) + ylim(0.5, 1) + theme(legend.position="none"),
+                  legend,
+                  align = 'vh',
+                  labels = c("a", "b", "c", "d", "e", "f", "g", "h"),
+                  hjust = -1,
+                  nrow = 3
+)
+
+pdf("/work-zfs/abattle4/parsana/networks_correction/publication_figures/suppfig8.pdf", height = 6.5, width = 7.2)
+print(fig2)
+dev.off()
+
diff --git a/publication_figures/suppfig15/average_module_size.png b/publication_figures/suppfig15/average_module_size.png
new file mode 100644
index 0000000..39fbd0c
Binary files /dev/null and b/publication_figures/suppfig15/average_module_size.png differ
diff --git a/publication_figures/suppfig15/gray_modules.R b/publication_figures/suppfig15/gray_modules.R
new file mode 100755
index 0000000..b19456a
--- /dev/null
+++ b/publication_figures/suppfig15/gray_modules.R
@@ -0,0 +1,44 @@
+library(reshape2)
+library(ggplot2)
+modsize = readRDS("wgcna_module_size.Rds")
+
+# hist(unlist(lapply(modsize[[1]], function(y) sapply(y, function(x) length(x)))), 
+#      breaks = 50, col = rgb(1, 0.75, 0.14, 0.5))
+# 
+# hist(unlist(lapply(modsize[[4]], function(y) length(y))), 
+#      col = rgb(0.39, 0.58, 0.92, 0.5), 
+#      add = T, breaks = 20)
+
+modraw <- modsize[[1]][c(46:50),]
+modpc <-  modsize[[4]][c(46:50),]
+
+modraw <- lapply(seq_len(ncol(modraw)), function(i) modraw[,i])
+modpc <- lapply(seq_len(ncol(modpc)), function(i) modpc[,i])
+
+grayraw <- sapply(modraw, function(x) sapply(x, function(y) y[1]))
+graypc <- sapply(modpc, function(x) sapply(x, function(y) y[1]))
+
+colnames(grayraw)<- colnames(modsize[[1]])
+colnames(graypc) <- colnames(modsize[[1]])
+
+grayraw <- melt(grayraw[,c(1:5)])
+graypc <- melt(graypc[,c(1:5)])
+
+grayraw$cutheights <- rep(seq(0.9,1.0,length.out = 50)[46:50], 5)
+grayraw$Type <- "uncorrected"
+graypc$cutheights <-  rep(seq(0.9,1.0,length.out = 50)[46:50], 5)
+graypc$Type <- "PC"
+graypc$Tissue <- graypc$Var2
+grayraw$Tissue <- grayraw$Var2
+
+graymerged <- rbind(grayraw, graypc)
+
+ggplot(graymerged, aes(x = cutheights, y = value, col = Type, shape = Tissue))+
+  geom_line(aes(linetype= Type))+
+  geom_point()+
+  xlab("Cut-heights")+
+  ylab("# of genes with unassigned module")+
+  theme_classic()
+
+ggsave(filename = "num_genes_gray_modules.png", width = 5, height = 4, units = "in", dpi = 600)
+
diff --git a/publication_figures/suppfig15/make_suppfig15.R b/publication_figures/suppfig15/make_suppfig15.R
new file mode 100755
index 0000000..c2171ce
--- /dev/null
+++ b/publication_figures/suppfig15/make_suppfig15.R
@@ -0,0 +1,9 @@
+library(cowplot)
+p1 <- ggdraw() + draw_image("average_module_size.png")
+p2 <- ggdraw() + draw_image("num_modules.png")
+p3 <- ggdraw() + draw_image("num_genes_gray_modules.png")
+p <-plot_grid(p1, p2, p3, ncol = 3, align = "h", labels = c("a", "b", "c"))
+
+pdf("../suppfig15.pdf", height = 3, width = 7)
+print(p)
+dev.off()
diff --git a/publication_figures/suppfig15/module_size.R b/publication_figures/suppfig15/module_size.R
new file mode 100755
index 0000000..d0cb66b
--- /dev/null
+++ b/publication_figures/suppfig15/module_size.R
@@ -0,0 +1,55 @@
+library(reshape2)
+library(ggplot2)
+modsize = readRDS("wgcna_module_size.Rds")
+
+# hist(unlist(lapply(modsize[[1]], function(y) sapply(y, function(x) length(x)))), 
+#      breaks = 50, col = rgb(1, 0.75, 0.14, 0.5))
+# 
+# hist(unlist(lapply(modsize[[4]], function(y) length(y))), 
+#      col = rgb(0.39, 0.58, 0.92, 0.5), 
+#      add = T, breaks = 20)
+
+modraw <- modsize[[1]][46:50,]
+modpc <-  modsize[[4]][46:50,]
+
+modraw <- lapply(seq_len(ncol(modraw)), function(i) modraw[,i])
+modpc <- lapply(seq_len(ncol(modpc)), function(i) modpc[,i])
+
+grayraw <- sapply(modraw, function(x) sapply(x, function(y) mean(y[-1])))
+graypc <- sapply(modpc, function(x) sapply(x, function(y) mean(y[-1])))
+colnames(grayraw)<- colnames(modsize[[1]])
+colnames(graypc) <- colnames(modsize[[1]])
+grayraw <- melt(grayraw[,c(1:5)])
+graypc <- melt(graypc[,c(1:5)])
+grayraw$cutheights <- rep(seq(0.9,1.0,length.out = 50)[46:50], 5)
+graypc$cutheights <-  rep(seq(0.9,1.0,length.out = 50)[46:50], 5)
+graypc$Tissue <- graypc$Var2
+grayraw$Tissue <- grayraw$Var2
+
+# grayraw <- sapply(modraw, function(x) unlist(lapply(x, function(y) y[-1])))
+# graypc <- sapply(modpc, function(x) unlist(lapply(x, function(y) y[-1])))
+# 
+# names(grayraw) <- colnames(modsize[[1]])
+# names(graypc) <- colnames(modsize[[4]])
+# 
+# grayraw <- melt(grayraw[1:5])
+# graypc <- melt(graypc[1:5])
+# 
+
+grayraw$Type <- "uncorrected"
+graypc$Type <- "PC"
+# graypc$Tissue <- graypc$L1
+# grayraw$Tissue <- grayraw$L1
+
+graymerged <- rbind(grayraw, graypc)
+
+ggplot(graymerged, aes(x = cutheights, y = value, col = Type,
+                       shape = Tissue))+
+  geom_line(aes(linetype= Type))+
+  geom_point()+
+  xlab("Cut-heights")+
+  ylab("average # genes per module")+
+  theme_classic()
+
+
+ggsave(filename = "average_module_size.png", width = 5, height = 4, units = "in", dpi = 600)
diff --git a/publication_figures/suppfig15/num_genes_gray_modules.png b/publication_figures/suppfig15/num_genes_gray_modules.png
new file mode 100644
index 0000000..f152f4a
Binary files /dev/null and b/publication_figures/suppfig15/num_genes_gray_modules.png differ
diff --git a/publication_figures/suppfig15/num_modules.R b/publication_figures/suppfig15/num_modules.R
new file mode 100755
index 0000000..b883edc
--- /dev/null
+++ b/publication_figures/suppfig15/num_modules.R
@@ -0,0 +1,45 @@
+library(reshape2)
+library(ggplot2)
+
+modsize = readRDS("wgcna_module_size.Rds")
+
+# hist(unlist(lapply(modsize[[1]], function(y) sapply(y, function(x) length(x)))), 
+#      breaks = 50, col = rgb(1, 0.75, 0.14, 0.5))
+# 
+# hist(unlist(lapply(modsize[[4]], function(y) length(y))), 
+#      col = rgb(0.39, 0.58, 0.92, 0.5), 
+#      add = T, breaks = 20)
+
+modraw <- modsize[[1]][c(46:50),]
+modpc <-  modsize[[4]][c(46:50),]
+
+modraw <- lapply(seq_len(ncol(modraw)), function(i) modraw[,i])
+modpc <- lapply(seq_len(ncol(modpc)), function(i) modpc[,i])
+
+grayraw <- sapply(modraw, function(x) sapply(x, function(y) length(y)-1))
+graypc <- sapply(modpc, function(x) sapply(x, function(y) length(y)-1))
+
+colnames(grayraw)<- colnames(modsize[[1]])
+colnames(graypc) <- colnames(modsize[[1]])
+ 
+grayraw <- melt(grayraw[,c(1:5)])
+graypc <- melt(graypc[,c(1:5)])
+
+grayraw$cutheights <- rep(seq(0.9,1.0,length.out = 50)[46:50], 5)
+grayraw$Type <- "uncorrected"
+graypc$cutheights <-  rep(seq(0.9,1.0,length.out = 50)[46:50], 5)
+graypc$Type <- "PC"
+graypc$Tissue <- graypc$Var2
+grayraw$Tissue <- grayraw$Var2
+
+graymerged <- rbind(grayraw, graypc)
+
+p <- ggplot(graymerged, aes(x = cutheights, y = value, col = Type, 
+                       shape = Tissue))+
+  geom_line(aes(linetype= Type))+
+  geom_point()+
+  xlab("Cut-heights")+
+  ylab("# of modules")+
+  theme_classic()
+
+ggsave(filename = "num_modules.png", plot = p, width = 5, height = 4, units = "in", dpi = 1500)
diff --git a/publication_figures/suppfig15/num_modules.png b/publication_figures/suppfig15/num_modules.png
new file mode 100644
index 0000000..abab9d2
Binary files /dev/null and b/publication_figures/suppfig15/num_modules.png differ
diff --git a/publication_figures/suppfig15/wgcna_module_size.Rds b/publication_figures/suppfig15/wgcna_module_size.Rds
new file mode 100755
index 0000000..e785224
Binary files /dev/null and b/publication_figures/suppfig15/wgcna_module_size.Rds differ
diff --git a/publication_figures/suppfig5/make_suppfig5.R b/publication_figures/suppfig5/make_suppfig5.R
new file mode 100644
index 0000000..2e99b41
--- /dev/null
+++ b/publication_figures/suppfig5/make_suppfig5.R
@@ -0,0 +1,80 @@
+library(parallel)
+library(recount)
+library(ggplot2)
+library(RColorBrewer)
+library(reshape2)
+library(cowplot)
+#### test association of expression PCs with gc coefficients
+
+  
+load("/work-zfs/abattle4/parsana/networks_correction/data/raw_subset.Rdata")
+load("/work-zfs/abattle4/parsana/networks_correction/data/pc_loadings.Rdata")
+
+  gc.coeff <- lapply(dat.expr, function(x) x$gc)
+  # fit linear model
+  gc.pc <- mapply(function(p,q,r){
+    lm.gc.fit <- lm(p[,1:r] ~ q)
+    lm.gc.fit
+    }, pc.loadings, gc.coeff, num.pc.estimates, SIMPLIFY = FALSE
+    )
+  # extract p-values
+  gc.pc.pvals <- lapply(gc.pc, function(lm.gc.fit){
+    pval <- sapply(summary(lm.gc.fit), function(x) x$coefficients[2,4])
+    pval <- p.adjust(pval, method = "BH")
+    pval
+    })
+  # extract percent variance explained/ R2
+  r2.pcs <- lapply(gc.pc, function(lm.gc.fit){
+    pval <- sapply(summary(lm.gc.fit), function(x) x$r.squared)
+    pval
+    })
+
+  # merge the lists to a matrix
+  r2.pcs <- sapply(r2.pcs, function(x) {length(x) <- 37; x})
+  gc.pc.pvals <- sapply(gc.pc.pvals, function(x) {length(x) <- 37; x})
+  rownames(r2.pcs) <- paste("PC",1:37,sep="")
+  rownames(gc.pc.pvals) <- paste("PC",1:37,sep="")
+
+# prepare data for plotting with ggplot2
+plot.r2 <- melt(r2.pcs)
+plot.pvals <- melt(gc.pc.pvals)
+
+  # set color pallete
+  # myPalette <- colorRampPalette(rev(brewer.pal(11, "Spectral")), space="Lab")
+  myPalette <- colorRampPalette(brewer.pal(9,'Blues'), space = "Lab")
+
+  # set default font size
+  theme_set(theme_classic(base_size = 9)) 
+  
+  # plot p-values
+  pval.fig <- ggplot(plot.pvals, aes(Var2, Var1)) + geom_tile(aes(fill = value),
+       colour = "white") + 
+    geom_text(size = 2, aes(label = prettyNum(value, digits=3, width=4, format="fg")
+)) +
+    scale_fill_gradientn(colours = rev(myPalette(10)[1:6]), limits = c(0,1)
+, values = c(0, 0.005, 0.05, 0.1, 0.5, 1)
+) + xlab("") + ylab("Principal Components") + theme(axis.text.x=element_text(colour="black", angle = 90), axis.text.y=element_text(colour="black"))
+  ggsave("pvals.png")
+
+  # plot r2
+  r2.fig <- ggplot(plot.r2, aes(Var2, Var1)) + geom_tile(aes(fill = value),
+       colour = "white") + 
+    geom_text(size = 2, aes(label = prettyNum(value, digits=3, width=4, format="fg")
+)) +
+    scale_fill_gradientn(colours = myPalette(10)[1:6], limits = c(0,1)
+    , values = c(0, 0.05, 0.1, 0.2, 0.3, 0.4, 1)
+) + xlab("") + ylab("Principal Components") + theme(axis.text.x=element_text(colour="black", angle = 90), axis.text.y=element_text(colour="black"))
+  ggsave("r2.png")
+
+fig.out <- plot_grid(pval.fig + theme(legend.position="bottom", legend.direction = "horizontal"),
+  r2.fig + theme(legend.position="bottom", legend.direction="horizontal"),
+  align = 'vh',
+           labels = c("a", "b"),
+           hjust = -1,
+           nrow = 1
+           )
+
+pdf("suppfig5.pdf", height = 8.5, width = 7)
+print(fig.out)
+dev.off()
+
diff --git a/publication_figures/suppfig5/suppfig5.pdf b/publication_figures/suppfig5/suppfig5.pdf
new file mode 100644
index 0000000..4bfee80
Binary files /dev/null and b/publication_figures/suppfig5/suppfig5.pdf differ