---
title: "Figure3"
author: "Liwen Xiao, Fengyi Zhang"
date: 2022-3
---

```{r load packages}
library(pheatmap)
library(WGCNA)
library(MASS)
library(ggplot2)
library(RColorBrewer)
library(reshape)
library(Rmisc)
library(MESS)
```



```{r set parameters}
rm(list = ls())

all.mod = 200
trans.m = 40
m1 = (all.mod - trans.m)/2
m2 = all.mod - m1
```



```{r generate simulated matrix}
add_cor <- function(eidim,max,ndim) {
  cor.vect <- c( )
  for (i in 1:ndim) {
    ei <- runif(eidim, -1, 1)
    cor.vect <- cbind(cor.vect, sqrt(max) * ei / sqrt(sum(ei^2) ) )
  }

  sole.cor <- t(cor.vect) %*% cor.vect
  diag(sole.cor) <- 1
  return(sole.cor)
}

sim.cor.runif.abs <- function(k,size,max,eidim) {
  ndim <- sum(size)
  bigcor<- matrix(rep(0, ndim*ndim), ncol=ndim)
  for (i in 1:(k) ){
    cor <- add_cor(eidim[i],max[i],size[i])
    if (i==1){bigcor[1:size[1], 1:size[1]] <- cor}
    if (i!=1){bigcor[(sum(size[1:(i-1)]) + 1):sum(size[1:i]),
                     (sum(size[1:(i-1)]) + 1):sum(size[1:i])] <- cor}
  }
  return(bigcor)
}
sim.cor <- abs(sim.cor.runif.abs(k = 2,size=c(m2,m1),max=c(0.5,0.5),eidim=c(2,2)))

add_noise <- function(eidim,epsilon,dim.row,dim.col) {
  bignoise <- matrix(nrow = dim.row,ncol = dim.col,0)
  for (j in 1:dim.col) {
    eivect <- c( )
    for (i in 1:(dim.row+1)) {
      ei <- runif(eidim, -1, 1)
      eivect <- cbind(eivect, sqrt(epsilon) * ei / sqrt(sum(ei^2) ) )
    }
    bigE <- t(eivect) %*% eivect
    noise_vect <- bigE[-1,1]
    bignoise[,j] <- noise_vect
  }
  return(bignoise)
}

noise_mat <- add_noise(eidim = 2,epsilon = 0.5,dim.row = m2,dim.col = m1)
sim.cor.noise <- sim.cor
sim.cor.noise[1:m2,(m2+1):all.mod] <- sim.cor.noise[1:m2,(m2+1):all.mod] - abs(noise_mat)
sim.cor.noise[(m2+1):all.mod,1:m2] <- sim.cor.noise[(m2+1):all.mod,1:m2] - t(abs(noise_mat))

###########################################################################################
rownames(sim.cor.noise) <- paste("taxon",seq(1,all.mod,1),sep="")
colnames(sim.cor.noise) <- paste("taxon",seq(1,all.mod,1),sep="")
adj.cor <- as.matrix(0.5+0.5*sim.cor.noise)

powers = seq(1,20,1)
library(igraph)
soft_thred <- function(mat,powers) {
  result <- data.frame(matrix(NA,length(powers),3))
  names(result) <- c("soft","r.squired","mean.connectivity")
  result$soft <- powers
  r.squired <- c()
  mean.connectivity <- c()
  for (pow in powers) {
    mat_soft <- mat
    for (i in 1:nrow(mat)) {
      for (j in 1:ncol(mat)) {
        mat_soft[i,j] <- mat[i,j]^pow
      }
    }
    mat_igraph <- graph_from_adjacency_matrix(mat_soft,mode = "undirected",weighted = T,diag = F)

    graph_density <- density(graph.strength(mat_igraph))
    y <- graph_density$y[intersect(which(graph_density$y != 0),which(graph_density$x != 0))]
    x <- graph_density$x[intersect(which(graph_density$y != 0),which(graph_density$x != 0))]
    r.squired <- c(r.squired,summary(lm(log(y)~log(x)))$r.squared)
    mean.connectivity <- c(mean.connectivity,sum(graph.strength(mat_igraph))/length(graph.strength(mat_igraph)))
  }
  result$r.squired <- r.squired
  result$mean.connectivity <- mean.connectivity
  return(result)
}

soft_cor <- soft_thred(adj.cor,powers)
pow_adjacency <- function(mat,pow) {
  mat_soft <- mat
  for (i in 1:nrow(mat)) {
    for (j in 1:ncol(mat)) {
      mat_soft[i,j] <- mat[i,j]^pow
    }
  }
  return(mat_soft)
}


soft_pow_cor <- 4
adj_soft_cor<- pow_adjacency(adj.cor,soft_pow_cor)
TOM_cor <- TOMsimilarity(adj_soft_cor)
dissTOM_cor <- 1-TOM_cor
# Call the hierarchical clustering function
hc_cor = hclust(as.dist(dissTOM_cor), method = "average");
# Plot the resulting clustering tree (dendrogram)
minModuleSize = m1
# Module identification using dynamic tree cut:
dynamicMods.cor = cutreeDynamic(dendro = hc_cor, distM = dissTOM_cor, deepSplit = 4, pamRespectsDendro = FALSE, minClusterSize = minModuleSize);
print(table(dynamicMods.cor))
dynamicColors.cor = labels2colors(dynamicMods.cor)
table(dynamicColors.cor) # Plot the dendrogram and colors underneath
####plotDendroAndColors(geneTree_zc, dynamicColors_zc, "Dynamic Tree Cut", dendroLabels = FALSE, hang = 0.03, addGuide = TRUE, guideHang = 0.05, main = "Gene dendrogram and module colors")
# Calculate topological overlap anew: this could be done more efficiently by saving the TOM
# calculated during module detection, but let us do it again here.
#dissTOM = 1-TOMsimilarityFromExpr(datExpr, power = 6);
# Transform dissTOM with a power to make moderately strong connections more visible in the heatmap
plotTOM_cor_num <- soft_pow_cor
plotTOM_cor = dissTOM_cor^plotTOM_cor_num;
# Set diagonal to NA for a nicer plot
diag(plotTOM_cor) = NA;
# Call the plot function
sim.cor.noise.p <- sim.cor.noise
diag(sim.cor.noise.p) <- 0
min_p <- min(sim.cor.noise.p)
max_p <- max(sim.cor.noise.p)
bk1 <- c(seq(-0.5,-0.05,by=0.05))
bk2 <- c(seq(0,0.5,by=0.05))
sim.cor.noise.p[sim.cor.noise.p > 0 & sim.cor.noise.p < 0.1] <- 0.1
sim.cor.noise.p[sim.cor.noise.p < 0 & sim.cor.noise.p > -0.1] <- -0.1
zc_row_mod <- data.frame(paste("m",dynamicMods.cor,sep=""))
rownames(zc_row_mod) <- rownames(sim.cor.noise.p)
names(zc_row_mod) <- "Mods"
zc_row_mod$Mods <- as.factor(zc_row_mod$Mods)
ann_colors = list(
  Mods = c(m1 = "#e34a33", m2 = "#2ca25f")
)



##plot fig3c heatmap
pheatmap(sim.cor.noise.p,cluster_rows = F,cluster_cols = F,
         color = c(colorRampPalette(colors = c("#034e7b","white"))(length(bk1)),colorRampPalette(colors = c("white","#990000"))(length(bk2))),
         breaks = c(bk1,bk2),show_rownames = F,show_colnames = F,annotation_row = zc_row_mod,annotation_col = zc_row_mod,annotation_colors = ann_colors[1],height=10,width=10)



```



```{r add random noise}
dist.cor.scale <- (dissTOM_cor - min(dissTOM_cor))/(max(dissTOM_cor)-min(dissTOM_cor))
new.cor <-as.matrix(dist.cor.scale)

distance_mds <- isoMDS(new.cor)
x = distance_mds$points[,1]
y = distance_mds$points[,2]
dat <- data.frame(x,y,dynamicMods.cor)
dat$dynamicMods.cor <- as.factor(dat$dynamicMods.cor)

add_descend_noise <- function(eidim,epsilon.max,epsilon.min,dim.row,dim.col){
  bignoise <- matrix(nrow = dim.row,ncol = dim.col,0)
  epsilon.step <- (epsilon.max-epsilon.min) / (dim.col-1)
  for (j in 1:dim.col) {
    epsilon <- epsilon.max - (j-1)*epsilon.step
    eivect <- c( )
    for (i in 1:(dim.row+1)) {
      ei <- runif(eidim, -1, 1)
      eivect <- cbind(eivect, sqrt(epsilon) * ei / sqrt(sum(ei^2) ) )
    }
    bigE <- t(eivect) %*% eivect
    noise_vect <- bigE[-1,1]
    bignoise[,j] <- noise_vect
  }
  return(bignoise)
}


descend_noise <- abs(add_descend_noise(eidim=2,epsilon.max = 1,epsilon.min = 0,dim.row = all.mod-m2, dim.col = trans.m))
range(descend_noise)
descend_noise_rev <- descend_noise
for (i in 1:ncol(descend_noise)) {
  descend_noise_rev[,i] <- descend_noise[,(ncol(descend_noise)-(i-1))]
}

#set.seed(123)
noise_order <- sample(1:trans.m,trans.m)
descend_noise_order <- descend_noise[,noise_order]
sim.cor.noise.descend <- sim.cor.noise
sim.cor.noise.descend[1:(m2-trans.m),(m2-trans.m+1):m2] <- sim.cor.noise.descend[1:(m2-trans.m),(m2-trans.m+1):m2] - descend_noise_order
sim.cor.noise.descend[(m2-trans.m+1):m2,1:(m2-trans.m)] <- sim.cor.noise.descend[(m2-trans.m+1):m2,1:(m2-trans.m)] - t(descend_noise_order)
sim.cor.noise.descend[(m2-trans.m+1):m2,(m2+1):all.mod] <- sim.cor.noise.descend[(m2-trans.m+1):m2,(m2+1):all.mod] + t(descend_noise_order)
sim.cor.noise.descend[(m2+1):all.mod,(m2-trans.m+1):m2] <- sim.cor.noise.descend[(m2+1):all.mod,(m2-trans.m+1):m2] + descend_noise_order

#####################################################################
adj.noise <- as.matrix(0.5+0.5*sim.cor.noise.descend)
soft_noise <- soft_thred(adj.noise,powers)
soft_pow_noise <- 4
adj_soft_noise<- pow_adjacency(adj.noise,soft_pow_noise)
TOM_noise <- TOMsimilarity(adj_soft_noise)
dissTOM_noise <- 1-TOM_noise
# Call the hierarchical clustering function
hc_noise = hclust(as.dist(dissTOM_noise), method = "average");
# Plot the resulting clustering tree (dendrogram)
minModuleSize = m1
# Module identification using dynamic tree cut:
dynamicMods.noise = cutreeDynamic(dendro = hc_noise, distM = dissTOM_noise, deepSplit = 4, pamRespectsDendro = FALSE, minClusterSize = minModuleSize);
print(table(dynamicMods.noise))
dynamicColors.noise = labels2colors(dynamicMods.noise)
table(dynamicColors.noise) # Plot the dendrogram and colors underneath
# Transform dissTOM with a power to make moderately strong connections more visible in the heatmap
plotTOM_noise_num <- soft_pow_noise
plotTOM_noise = dissTOM_noise^plotTOM_noise_num;
# Set diagonal to NA for a nicer plot
diag(plotTOM_noise) = NA;
# Call the plot function
if (dynamicMods.noise[1] == 2) {
  dynamicMods.noise[which(dynamicMods.noise == 1)] <- 0
  dynamicMods.noise[which(dynamicMods.noise == 2)] <- 1
  dynamicMods.noise[which(dynamicMods.noise == 0)] <- 2
}

sim.cor.noise.descend.p <- sim.cor.noise.descend
diag(sim.cor.noise.descend.p) <- 0
min_p <- min(sim.cor.noise.descend.p)

max_p <- max(sim.cor.noise.descend.p)


bk1 <- c(seq(-0.5,-0.05,by=0.05))
bk2 <- c(seq(0,0.5,by=0.05))
sim.cor.noise.descend.p[sim.cor.noise.descend.p > 0 & sim.cor.noise.descend.p < 0.1] <- 0.1
sim.cor.noise.descend.p[sim.cor.noise.descend.p < 0 & sim.cor.noise.descend.p > -0.1] <- -0.1

sim.cor.noise.descend.p[sim.cor.noise.descend.p > 0.5] <- 0.5
sim.cor.noise.descend.p[sim.cor.noise.descend.p < -0.5] <- -0.5
zc_row_mod <- data.frame(paste("m",dynamicMods.noise,sep=""))
rownames(zc_row_mod) <- rownames(sim.cor.noise.descend.p)
names(zc_row_mod) <- "Mods"
zc_row_mod$Mods <- as.factor(zc_row_mod$Mods)
ann_colors = list(
  Mods = c(m1 = "#e34a33", m2 = "#2ca25f")
)


###plot fig3c disturbed matrix
pheatmap(sim.cor.noise.descend.p,cluster_rows = F,cluster_cols = F,
         color = c(colorRampPalette(colors = c("#034e7b","white"))(length(bk1)),colorRampPalette(colors = c("white","#990000"))(length(bk2))),
         breaks = c(bk1,bk2),show_rownames = F,show_colnames = F,annotation_row = zc_row_mod,annotation_col = zc_row_mod,annotation_colors = ann_colors[1],height=10,width=10)

as.numeric(dynamicMods.cor)
as.numeric(dynamicMods.noise)
dynamicMods_inter <- paste("control_mod_",dynamicMods.cor,"-","case_mod_",dynamicMods.noise,sep="")
dynamicMods_inter_names <- names(table(dynamicMods_inter))

dist.noise.scale <- (dissTOM_noise - min(dissTOM_noise))/(max(dissTOM_noise)-min(dissTOM_noise))

```



```{r Figure 3e}
mynode = read.table("node_result_all.txt",header = T,sep = '\t')
mynode$jaccard.scale = 1-mynode$jaccard
mynode$nesh.scale = (mynode$nesh-min(mynode$nesh))/(max(mynode$nesh)-min(mynode$nesh))
mynode$netzgo.scale = (mynode$netzgo-min(mynode$netzgo))/(max(mynode$netzgo)-min(mynode$netzgo))

#table(mynode[,c(5,6)])
mod = c(60,60,40,60,80,80,50,50,40,70)
for (k in 1:10)
{
  a = mynode[which(mynode$perm == (k + 10)),]
  rownames(a) = a$taxon_names
  #b = melt(a[,c(1,7,8,9)],id.vars = "taxon_names")
  b = melt(a[,c(1,2,3,4)],id.vars = "taxon_names")
  b$class = a[as.character(b$taxon_names),"class"]
  b$tax = substr(b$taxon_names,6,8)
  b$variable = as.character(b$variable)
  b[which(b$variable == "netzgo"),"variable"] = "NetMoss"
  b[which(b$variable == "nesh"),"variable"] = "Nesh"
  b[which(b$variable == "jaccard"),"variable"] = "Jaccard"
  b$variable = factor(b$variable,levels = c("NetMoss","Nesh","Jaccard"))
  #density
  p1 = ggplot(b,aes(x = value,fill = class))+
    geom_density(alpha = .5)+
    facet_grid(.~variable,scales = "free")+
    #scale_fill_brewer(palette = "Set1")+
    scale_fill_manual(values = c("#989898","#CD524B"))+
    theme_bw()+
    theme(panel.grid = element_blank())
  #point
  b$tax = factor(b$tax,levels = b[1:200,5])
  p2 = ggplot(b,aes(tax,value,color = class))+
    geom_point(size = .5)+coord_flip()+
    facet_grid(.~variable,scales = "free")+
    #scale_color_brewer(palette = "Set1")+
    scale_color_manual(values = c("#989898","#CD524B"))+
    theme_bw()+
    theme(panel.grid = element_blank())+
    theme(axis.text.y = element_blank(),axis.ticks.y = element_blank())+
    labs(x = "Taxon")+
    geom_rect(aes(xmin=((200-mod[k])/2), xmax=((200+mod[k])/2+1), ymin=-Inf, ymax=Inf),fill = "grey60",alpha = .01, color = NA)+
    geom_point(size = .5)
  eval(parse(text = paste0("p",k,"1 = ggplot(b,aes(x = value,fill = class))+
                              geom_density(alpha = .5)+
                              facet_grid(.~variable,scales = 'free')+
                              scale_fill_manual(values = c('#989898','#CD524B'),guide=FALSE)+
                              theme_bw()+
                              labs(x = 'Score', y = 'Density')+
                              theme(panel.grid = element_blank(),)")))
  ##point all
  eval(parse(text = paste0("p",k,"2 = ggplot(b,aes(tax,value,color = class))+
                             geom_point(size = .05)+coord_flip()+
                             facet_grid(.~variable,scales = 'free')+
                             scale_color_manual(values = c('#989898','#CD524B'),guide=FALSE)+
                             theme_bw()+
                             theme(panel.grid = element_blank())+
                             theme(axis.text.y = element_blank(),axis.ticks.y = element_blank())+
                             labs(x = 'Taxon',y = 'Score')+
                             geom_rect(aes(xmin=((200-mod[",k,"])/2), xmax=((200+mod[",k,"])/2), ymin=-Inf, ymax=Inf),fill = 'grey60',alpha = .01, color = NA)+
                             geom_point(size = .5)")))
}


#net threshold
net = b[which(b$variable == "NetMoss"),]
net.target = net[which(net$class == "yes"),]
y.net = data.frame()
x = 0
for (k in 1:100)
{
  y.net[k,1] = x
  y.net[k,2] = nrow(net.target[which(net.target$value > x),])
  y.net[k,3] = nrow(net[which(net$value > x),])-nrow(net.target[which(net.target$value > x),])
  y.net[k,4] = nrow(net.target)-nrow(net.target[which(net.target$value > x),])
  y.net[k,5] = nrow(net)-nrow(net[which(net$value > x),])-nrow(net.target)+nrow(net.target[which(net.target$value > x),])
  y.net[k,6] = y.net[k,2]/(y.net[k,2]+y.net[k,4])
  y.net[k,7] = y.net[k,5]/(y.net[k,3]+y.net[k,5])
  x = x + max(net$value)/100
}
colnames(y.net) = c("threshold","TP","FP","FN","TN","sensitivity","specificity")
y2.net = melt(y.net[,c(1,6,7)],id.vars = "threshold")

#nesh threshold
nesh = b[which(b$variable == "Nesh"),]
nesh.target = nesh[which(nesh$class == "yes"),]
y.nesh = data.frame()
x = 0
for (k in 1:100)
{
  y.nesh[k,1] = x
  y.nesh[k,2] = nrow(nesh.target[which(nesh.target$value > x),])
  y.nesh[k,3] = nrow(nesh[which(nesh$value > x),])-nrow(nesh.target[which(nesh.target$value > x),])
  y.nesh[k,4] = nrow(nesh.target)-nrow(nesh.target[which(nesh.target$value > x),])
  y.nesh[k,5] = nrow(nesh)-nrow(nesh[which(nesh$value > x),])-nrow(nesh.target)+nrow(nesh.target[which(nesh.target$value > x),])
  y.nesh[k,6] = y.nesh[k,2]/(y.nesh[k,2]+y.nesh[k,4])
  y.nesh[k,7] = y.nesh[k,5]/(y.nesh[k,3]+y.nesh[k,5])
  x = x + max(nesh$value)/100
}
colnames(y.nesh) = c("threshold","TP","FP","FN","TN","sensitivity","specificity")
y2.nesh = melt(y.nesh[,c(1,6,7)],id.vars = "threshold")

#jaccard threshold
jaccard = b[which(b$variable == "Jaccard"),]
jaccard.target = jaccard[which(jaccard$class == "yes"),]
y.jaccard = data.frame()
x = 0
for (k in 1:100)
{
  y.jaccard[k,1] = x
  y.jaccard[k,2] = nrow(jaccard.target[which(jaccard.target$value < x),])
  y.jaccard[k,3] = nrow(jaccard[which(jaccard$value < x),])-nrow(jaccard.target[which(jaccard.target$value < x),])
  y.jaccard[k,4] = nrow(jaccard.target)-nrow(jaccard.target[which(jaccard.target$value < x),])
  y.jaccard[k,5] = nrow(jaccard)-nrow(jaccard[which(jaccard$value < x),])-nrow(jaccard.target)+nrow(jaccard.target[which(jaccard.target$value < x),])
  y.jaccard[k,6] = y.jaccard[k,5]/(y.jaccard[k,3]+y.jaccard[k,5])
  y.jaccard[k,7] = y.jaccard[k,2]/(y.jaccard[k,2]+y.jaccard[k,4])
  x = x + max(jaccard$value)/100
}
colnames(y.jaccard) = c("threshold","TP","FP","FN","TN","sensitivity","specificity")
y2.jaccard = melt(y.jaccard[,c(1,6,7)],id.vars = "threshold")

#roc
y.net.all = data.frame()
y.nesh.all = data.frame()
y.jaccard.all = data.frame()
mod = c(60,60,40,60,80,80,50,50,40,70)
for (i in 1:10)
{
  a = mynode[which(mynode$perm == (i + 10)),]
  rownames(a) = a$taxon_names
  b = melt(a[,c(1,2,3,4)],id.vars = "taxon_names")
  b$class = a[as.character(b$taxon_names),"class"]
  b$tax = substr(b$taxon_names,6,8)
  d = b[which(as.numeric(b$tax) > ((200-mod[i])/2)),]
  d = d[which(as.numeric(d$tax) < ((200+mod[i])/2+1)),]
  #net threshold
  net = b[which(b$variable == "netzgo"),]
  net.target = net[which(net$class == "yes"),]
  y.net = data.frame()
  x = 0
  for (k in 1:100)
  {
    y.net[k,1] = x
    y.net[k,2] = nrow(net.target[which(net.target$value > x),])
    y.net[k,3] = nrow(net[which(net$value > x),])-nrow(net.target[which(net.target$value > x),])
    y.net[k,4] = nrow(net.target)-nrow(net.target[which(net.target$value > x),])
    y.net[k,5] = nrow(net)-nrow(net[which(net$value > x),])-nrow(net.target)+nrow(net.target[which(net.target$value > x),])
    y.net[k,6] = y.net[k,2]/(y.net[k,2]+y.net[k,4])
    y.net[k,7] = y.net[k,5]/(y.net[k,3]+y.net[k,5])
    x = x + max(net$value)/100
  }
  colnames(y.net) = c("threshold","TP","FP","FN","TN","sensitivity","specificity")
  #nesh threshold
  nesh = d[which(d$variable == "nesh"),]
  nesh.target = nesh[which(nesh$class == "yes"),]
  y.nesh = data.frame()
  x = 0
  for (k in 1:100)
  {
    y.nesh[k,1] = x
    y.nesh[k,2] = nrow(nesh.target[which(nesh.target$value > x),])
    y.nesh[k,3] = nrow(nesh[which(nesh$value > x),])-nrow(nesh.target[which(nesh.target$value > x),])
    y.nesh[k,4] = nrow(nesh.target)-nrow(nesh.target[which(nesh.target$value > x),])
    y.nesh[k,5] = nrow(nesh)-nrow(nesh[which(nesh$value > x),])-nrow(nesh.target)+nrow(nesh.target[which(nesh.target$value > x),])
    y.nesh[k,6] = y.nesh[k,2]/(y.nesh[k,2]+y.nesh[k,4])
    y.nesh[k,7] = y.nesh[k,5]/(y.nesh[k,3]+y.nesh[k,5])
    x = x + max(nesh$value)/100
  }
  colnames(y.nesh) = c("threshold","TP","FP","FN","TN","sensitivity","specificity")
  #jaccard threshold
  jaccard = d[which(d$variable == "jaccard"),]
  jaccard.target = jaccard[which(jaccard$class == "yes"),]
  y.jaccard = data.frame()
  x = 0
  for (k in 1:100)
  {
    y.jaccard[k,1] = x
    y.jaccard[k,2] = nrow(jaccard.target[which(jaccard.target$value < x),])
    y.jaccard[k,3] = nrow(jaccard[which(jaccard$value < x),])-nrow(jaccard.target[which(jaccard.target$value < x),])
    y.jaccard[k,4] = nrow(jaccard.target)-nrow(jaccard.target[which(jaccard.target$value < x),])
    y.jaccard[k,5] = nrow(jaccard)-nrow(jaccard[which(jaccard$value < x),])-nrow(jaccard.target)+nrow(jaccard.target[which(jaccard.target$value < x),])
    y.jaccard[k,6] = y.jaccard[k,2]/(y.jaccard[k,2]+y.jaccard[k,4])
    y.jaccard[k,7] = y.jaccard[k,5]/(y.jaccard[k,3]+y.jaccard[k,5])
    x = x + max(jaccard$value)/100
  }
  colnames(y.jaccard) = c("threshold","TP","FP","FN","TN","sensitivity","specificity")
  y.net$auc = auc(1-y.net$specificity,y.net$sensitivity,type = 'spline')
  y.nesh$auc = auc(1-y.nesh$specificity,y.nesh$sensitivity,type = 'spline')
  y.jaccard$auc = auc(1-y.jaccard$specificity,y.jaccard$sensitivity,type = 'spline')
  y.net$perm = (i+10)
  y.nesh$perm = (i+10)
  y.jaccard$perm = (i+10)
  y.net.all = rbind(y.net,y.net.all)
  y.nesh.all = rbind(y.nesh,y.nesh.all)
  y.jaccard.all = rbind(y.jaccard,y.jaccard.all)
}
y.net.all$type = "NetMoss"
y.nesh.all$type = "Nesh"
y.jaccard.all$type = "Jaccard"
y.all = rbind(y.nesh.all,y.net.all,y.jaccard.all)

####################plot fig3e roc###########################
ggplot(y.all,aes(1-specificity,sensitivity,color = factor(perm)))+
  geom_abline(slope = 1)+
  geom_path(size = 1)+
  scale_y_continuous(limits = c(0,1))+
  facet_grid(.~type)+
  scale_color_manual(guide = "none",values = colorRampPalette(brewer.pal( 8, "Set2"))(10))+
  theme_bw()+
  theme(panel.grid = element_blank())+
  labs(x = "False positive frcation", y = "True positive fraction")




```

