---
title: "Figure1"
author: "Liwen Xiao, Fengyi Zhang"
date: "2022-3"
output: pdf_document
---

```{r load packages}
library(tidyverse)
library(coin)
library(pROC)
library(yaml)
library(RColorBrewer)
library(cowplot)
library(vegan)
library(labdsv)
library(plotrix)
library(ggplot2)
library(igraph)
```



```{r read crc files}
rm(list = ls())
crc_num <- c(1,2,3,4,5,6,7)
crc_feat_list <- list()
zc_feat_list <- list()
for (i in crc_num) {
  #i <- 1
  crc_feat <- read.table(paste0("d_crc",i,".txt"),header = T,sep="\t")
  zc_feat <- read.table(paste0("h_crc",i,".txt"),header = T,sep="\t")
  rownames(crc_feat) <- as.character(crc_feat$X)
  rownames(zc_feat) <- as.character(zc_feat$X)
  rownames(crc_feat) <- gsub("\\.","-",rownames(crc_feat))
  colnames(crc_feat) <- gsub("\\.","-",colnames(crc_feat))
  rownames(crc_feat) <- gsub(" ","-",rownames(crc_feat))
  colnames(crc_feat) <- gsub(" ","-",colnames(crc_feat))
  rownames(crc_feat) <- gsub('\\[','',rownames(crc_feat))
  colnames(crc_feat) <- gsub('\\[','',colnames(crc_feat))
  rownames(crc_feat) <- gsub('\\]','',rownames(crc_feat))
  colnames(crc_feat) <- gsub('\\]','',colnames(crc_feat))
  rownames(zc_feat) <- gsub("\\.","-",rownames(zc_feat))
  colnames(zc_feat) <- gsub("\\.","-",colnames(zc_feat))
  rownames(zc_feat) <- gsub(" ","-",rownames(zc_feat))
  colnames(zc_feat) <- gsub(" ","-",colnames(zc_feat))
  rownames(zc_feat) <- gsub('\\[','',rownames(zc_feat))
  colnames(zc_feat) <- gsub('\\[','',colnames(zc_feat))
  rownames(zc_feat) <- gsub('\\]','',rownames(zc_feat))
  colnames(zc_feat) <- gsub('\\]','',colnames(zc_feat))
  crc_feat <- crc_feat[,-1]
  zc_feat <- zc_feat[,-1]
  crc_feat_list[[i]] <- crc_feat
  zc_feat_list[[i]] <- zc_feat
}

crc_feat_all <- as.data.frame(crc_feat_list[[crc_num[1]]])
crc_feat_all$id <- rownames(crc_feat_all)
for (i in crc_num[-1]) {
  try <- as.data.frame(crc_feat_list[[i]])
  try$id <- rownames(try)
  crc_feat_all <- merge(crc_feat_all,try,by.x ="id",all=T)
}

zc_feat_all <- as.data.frame(zc_feat_list[[crc_num[1]]])
zc_feat_all$id <- rownames(zc_feat_all)
for (i in crc_num[-1]) {
  try <- as.data.frame(zc_feat_list[[i]])
  try$id <- rownames(try)
  zc_feat_all <- merge(zc_feat_all,try,by.x ="id",all=T)
}
g.int = intersect(crc_feat_all$id,zc_feat_all$id)
crc_feat_all = crc_feat_all[as.character(crc_feat_all$id) %in% as.character(g.int),]
zc_feat_all = zc_feat_all[as.character(zc_feat_all$id) %in% as.character(g.int),]
zc_feat_all <- zc_feat_all[,-1]
crc_feat_all <- cbind(crc_feat_all,zc_feat_all)
crc_feat_all[is.na(crc_feat_all)] <- 0
rownames(crc_feat_all) <- crc_feat_all$id
crc_feat_all <- crc_feat_all[,-1]
crc_feat_all <- as.matrix(crc_feat_all)


union_genus <- rownames(crc_feat_list[[1]])
for (i in crc_num) {
  eval(parse(text = paste("union_genus <- union(union_genus,rownames(crc_feat_list[[",i,"]]))",sep="")))
  
}

union_genus2 <- rownames(zc_feat_list[[1]])
for (i in crc_num) {
  eval(parse(text = paste("union_genus2 <- union(union_genus2,rownames(zc_feat_list[[",i,"]]))",sep="")))
  
}

union_genus = unique(c(union_genus,union_genus2))

for (i in crc_num) {
  crc_feat <- crc_feat_list[[i]]
  crc_feat_union <- matrix(NA, nrow=length(union_genus), ncol=ncol(crc_feat), 
                           dimnames=list(union_genus, colnames(crc_feat)))
  crc_feat_union[rownames(crc_feat),] <- as.matrix(crc_feat)
  crc_feat_union[is.na(crc_feat_union)] <- 0
  zc_feat <- zc_feat_list[[i]]
  zc_feat_union <- matrix(NA, nrow=length(union_genus), ncol=ncol(zc_feat), 
                          dimnames=list(union_genus, colnames(zc_feat)))
  zc_feat_union[rownames(zc_feat),] <- as.matrix(zc_feat)
  zc_feat_union[is.na(zc_feat_union)] <- 0
  crc_feat_list[[i]] <- crc_feat_union
  zc_feat_list[[i]] <- zc_feat_union
}

p.val <- matrix(NA, nrow=length(union_genus), ncol=(length(crc_num)+1), 
                dimnames=list(union_genus, c(paste0("crc",crc_num), 'all')))
fc <- p.val

tag <- "species"
parameters <- yaml.load_file(paste0("parameters.yaml"))
log.n0 <- ifelse(tag=='species',
                 as.numeric(parameters$associations$log.n0),
                 as.numeric(parameters$associations$log.n0.func))
mult.corr <- parameters$associations$mult.corr

cat("Calculating effect size for every feature...\n")
pb <- txtProgressBar(max=length(union_genus), style=3)

# caluclate wilcoxon test and effect size for each feature and study
for (f in union_genus) {
  # for each study
  for (s in crc_num) {
    x <- crc_feat_list[[s]][f,]
    y <- zc_feat_list[[s]][f,]
    which(x!=0)
    which(y!=0)
    #wilcox.test(x, y, exact=FALSE)$p.value
    # Wilcoxon
    p.val[f,paste0("crc",s)] <- wilcox.test(x, y, exact=FALSE)$p.value
    
    # FC
    q.p <- quantile(log10(x+log.n0), probs=seq(.1, .9, .05))
    q.n <- quantile(log10(y+log.n0), probs=seq(.1, .9, .05))
    fc[f, paste0("crc",s)] <- sum(q.p - q.n)/length(q.p)
  }
  y <- c()
  x <- c()
  block <- c()
  for (s in crc_num) {
    y <- c(y,crc_feat_list[[s]][f,],zc_feat_list[[s]][f,])
    x <- c(x,rep("crc",length(crc_feat_list[[s]][f,])),rep("zc",length(zc_feat_list[[s]][f,])))
    block <- c(block,rep(paste0("crc",s),(length(crc_feat_list[[s]][f,])+length(zc_feat_list[[s]][f,]))))
  }
  d <- data.frame(y = y,x=x,block=block)
  p.val[f,'all'] <- pvalue(wilcox_test(y ~ x | block, data=d))
  
  # FC
  fc[f, 'all'] <- mean(fc[f, c(paste0("crc",crc_num))])
  
  # progressbar
  setTxtProgressBar(pb, (pb$getVal()+1))
}
cat('\n')

p.adj <- data.frame(apply(p.val, MARGIN=2, FUN=p.adjust, method=mult.corr),
                    check.names = FALSE)

```



```{r Figure 1b} 
crc_feat_all <- crc_feat_all[,-which(colSums(crc_feat_all) == 0)]

dist = vegdist(t(crc_feat_all), method = 'bray')
pco.results = pco(dist, k=2)

axis.1.title <- paste('PCo1 [', 
                      round((pco.results$eig[1]/sum(pco.results$eig))*100,1),
                      '%]', sep='')
axis.2.title <- paste('PCo2 [', 
                      round((pco.results$eig[2]/sum(pco.results$eig))*100,1),
                      '%]', sep='')

meta_all <- read.csv(paste("ZC.csv",sep = ""),as.is=T) 
meta_crc <- meta_all[meta_all$sample.ID %in% rownames(pco.results$points),]
table(meta_crc$dataset)
table(meta_crc$gender)

meta_crc_order <- c()
for (i in rownames(pco.results$points)){
  meta_crc_order <- c(meta_crc_order,which(meta_crc$sample.ID == i))
}
meta_crc <- meta_crc[meta_crc_order,]


df.plot <- tibble(Axis1 = -1*pco.results$points[,1],
                  Axis2 = pco.results$points[,2],
                  Sample_ID = rownames(pco.results$points),
                  Group=meta_crc$gender,
                  Study=meta_crc$dataset)

study.cols <- c("#2ca25f","#8856a7","#e34a33","#3182bd","#F2CC30","#737373","#78c679")
names(study.cols) <- paste0("crc",crc_num)
group.cols <- c("#C2C2C2","#d9544d")
names(group.cols) <- c("zc","de")
# ##############################################################################
# subplots

# main plot 
g.main <- df.plot %>% 
  ggplot(aes(x=Axis1, y=Axis2, shape=Group, col=Study)) +
  geom_point() + 
  scale_colour_manual(values=study.cols, guide=FALSE) + 
  scale_shape_manual(values=c(19, 1), guide=FALSE) + 
  scale_x_continuous(position='top') +
  xlab(axis.1.title) + ylab(axis.2.title) +
  theme(panel.background = element_rect(fill='white', color = 'black'),
        axis.ticks=element_blank(), axis.text = element_blank(),
        axis.title.x = element_blank(),axis.title.y = element_blank(),
        panel.grid = element_blank())

#g.main
# study boxplot axis 1
g.s.1 <- df.plot %>% 
  mutate(Study=factor(Study, levels=names(study.cols))) %>% 
  ggplot(aes(y=Axis1, x=Study, fill=Study)) + 
  geom_boxplot(size=0.1) + 
  scale_fill_manual(values=study.cols, guide=FALSE) +
  theme(axis.ticks = element_blank(),
        panel.background = element_rect(fill='white', color = 'black'),
        axis.text = element_blank(), 
        axis.title.x = element_blank(),axis.title.y = element_blank(),
        panel.grid = element_blank()) + 
  coord_flip()
# study boxplot axis 2
g.s.2 <- df.plot %>% 
  mutate(Study=factor(Study, levels=names(study.cols))) %>% 
  ggplot(aes(y=Axis2, x=Study, fill=Study)) + 
  geom_boxplot(size=0.1) + 
  scale_fill_manual(values=study.cols, guide=FALSE) +
  scale_x_discrete(position='top') +
  theme(axis.ticks=element_blank(), 
        panel.background = element_rect(fill='white', color = 'black'),
        axis.text = element_blank(),
        axis.title.y = element_blank(),
        axis.title.x = element_blank(),
        panel.grid = element_blank())
# group plot axis1
g.g.1 <- df.plot %>% 
  ggplot(aes(x=Group, y=Axis1, fill=Group)) +
  geom_boxplot(size=0.1) +
  scale_fill_manual(values=group.cols, guide=FALSE) + 
  ylab(axis.1.title) + 
  theme(axis.ticks.y=element_blank(),
        axis.text.y=element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_blank(),
        panel.background = element_rect(fill='white', color='black'),
        panel.grid = element_blank()) + 
  coord_flip()

# group plot axis2
g.g.2 <- df.plot %>% 
  ggplot(aes(x=Group, y=Axis2, fill=Group)) +
  geom_boxplot(size=0.1) +
  scale_fill_manual(values=group.cols, guide=FALSE) + 
  scale_x_discrete(position='top') + 
  scale_y_continuous(position = 'right') +
  ylab(axis.2.title) + 
  theme(axis.ticks.x=element_blank(),
        axis.text.x=element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        panel.background = element_rect(fill='white', color='black'),
        panel.grid = element_blank())

# ##############################################################################
### Plot  fig 1b
#pdf("fig1b.pdf",width = 5,height = 5)
plot_grid(g.main, g.s.2, g.g.2, g.s.1, NULL, NULL, g.g.1, NULL, NULL,
          nrow=3,
          rel_widths = c(0.8, 0.2, 0.2), rel_heights = c(0.8, 0.2, 0.2))
#dev.off()
```



```{r Figure 1d}
theme_set(theme_gray())
p.vals <- as.data.frame(p.adj)
fc.mat <- as.data.frame(fc)
ref.studies <- paste0("crc",crc_num)

#species.heatmap <- 
#  rownames(p.vals)[which(p.vals$all < 
#                           parameters$associations$alpha.single.study)]

species.heatmap <-
  rownames(p.vals)[which(p.vals$all < 0.05)]

fc.sign <- sign(fc.mat)
fc.sign[fc.sign == 0] <- 1

p.val.signed <- -log10(p.vals[species.heatmap,"all", drop=FALSE])*
  fc.sign[species.heatmap, 'all']
top.markers <- rownames(p.val.signed[is.infinite(p.val.signed$all),,drop=FALSE])
p.val.signed[top.markers, 'all'] <- 100 

#species.heatmap.orderd <- rownames(p.val.signed[order(p.val.signed$all),,
#                                                drop=FALSE])

p.val.signed$q = p.vals[species.heatmap,"all"]
species.heatmap.orderd <- rownames(p.val.signed[order(p.val.signed$q,decreasing = T),,drop=FALSE])

# take only those
fc.mat.plot <- fc.mat[species.heatmap.orderd,]
p.vals.plot <- p.vals[species.heatmap.orderd,]

# ##############################################################################
# prepare plotting

# colorscheme for fc heatmap 
mx <- max(abs(range(fc.mat.plot, na.rm=TRUE)))
mx <- ifelse(round(mx, digits = 1) < mx, 
             round(mx, digits = 1) + 0.1, 
             round(mx, digits = 1))
brs = seq(-mx, mx, by=0.05)
num.col.steps = length(brs)-1
n = floor(0.45*num.col.steps)
col.hm = c(rev(colorRampPalette(brewer.pal(9, 'Blues'))(n)),
           rep('#FFFFFF', num.col.steps-2*n),
           colorRampPalette(brewer.pal(9, 'Reds'))(n))
# color scheme for pval heatmap
alpha.breaks=c(1e-06, 1e-05, 1e-04, 1e-03, 1e-02, 1e-01)
p.vals.bin <- data.frame(apply(p.vals.plot, 2, FUN=.bincode, 
                               breaks = c(0, alpha.breaks, 1), 
                               include.lowest = TRUE),
                         check.names = FALSE)
p.val.greys <-  c(paste0('grey', 
                         round(seq(from=10, to=80, 
                                   length.out = length(alpha.breaks)))), 
                  'white')
names(p.val.greys) <- as.character(1:7)

bar_plot <- data.frame(X = 1:length(brs),Y = brs)
bar_plot$X <- as.factor(bar_plot$X)
g1 <- ggplot(bar_plot,aes(x=X, y=1, fill=Y)) + 
  geom_tile() + theme_minimal() + 
  theme(axis.text = element_blank(),
        axis.ticks = element_blank(),
        axis.title = element_blank(), 
        panel.grid = element_blank(),
        panel.background = element_rect(fill=NULL, colour='black'),
        plot.margin = unit(c(0, 0, 0, 0), 'cm')) + 
  scale_y_continuous(expand = c(0, 0)) + 
  scale_fill_gradientn(colours=col.hm, limits=c(-mx, mx), guide=FALSE)
# function to plot both into a grid
plot.single.study.heatmap <- function(x){
  df.plot <- tibble(species=factor(rownames(p.vals.plot), 
                                   levels=rev(rownames(p.vals.plot))),
                    p.vals=as.factor(p.vals.bin[[x]]),
                    fc=fc.mat.plot[[x]])
  g1 <- df.plot %>% 
    ggplot(aes(x=species, y=1, fill=fc)) + 
    geom_tile() + theme_minimal() + 
    theme(axis.text = element_blank(),
          axis.ticks = element_blank(),
          axis.title = element_blank(), 
          panel.grid = element_blank(),
          panel.background = element_rect(fill=NULL, colour='black'),
          plot.margin = unit(c(0, 0, 0, 0), 'cm')) + 
    scale_y_continuous(expand = c(0, 0)) + 
    scale_fill_gradientn(colours=col.hm, limits=c(-mx, mx), guide=FALSE)
}

# ##############################################################################
# plot
# p.value histogram
g1 <- tibble(species=factor(rownames(p.vals.plot), 
                            levels=rev(rownames(p.vals.plot))),
             p.vals=-log10(p.vals.plot$all),
             colour=p.vals > 5) %>% 
  ggplot(aes(x=species, y=p.vals, fill=colour)) + 
  geom_bar(stat='identity') + 
  theme_classic() + 
  xlab('Gut microbial genera') + 
  ylab('-log10(q-value)') + 
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks.x = element_blank(),
        axis.text.x = element_text(angle = 90),
        panel.background = element_rect(fill=NULL, colour = 'black')) + 
  scale_y_continuous(limits=c(0, 15), expand = c(0, 0)) + 
  scale_x_discrete(position='top') + 
  scale_fill_manual(values=c('lightgrey', 'darkgrey'), guide=FALSE) +
  geom_hline(aes(yintercept = 5),colour="#990000", linetype="dashed")

g.lst <- lapply(ref.studies, plot.single.study.heatmap)

##plot fig1d
plot_grid(g1, g.lst[[1]], g.lst[[2]], g.lst[[3]], g.lst[[4]],g.lst[[5]],g.lst[[6]],g.lst[[7]],
          ncol=1, align = 'v', rel_heights = c(0.44,rep(0.08, 7)))

```


```{r Figure 1c}
colnames(p.adj)
for (i in crc_num) {
  num <- which(crc_num == i)
  print(length(which(p.adj[,num] < 0.05)))
  
}

for (i in crc_num) {
  num <- which(crc_num == i)
  print(length(which(p.val[,num] < 0.05)))
  
}

pval_samples <- c()
pval_genus_list <- list()
padj_samples <- c()
padj_genus_list <- list()
for (i in crc_num) {
  num <- which(crc_num == i)
  crc_genus <- union_genus[which(p.val[,num] < 0.05)]
  pval_samples  <- c(pval_samples,length(crc_genus))
  pval_genus_list[[i]] <-  crc_genus
  crc_genus <- union_genus[which(p.adj[,num] < 0.05)]
  padj_samples  <- c(padj_samples,length(crc_genus))
  padj_genus_list[[i]] <-  crc_genus
  
}

crc_pval_genus <- pval_genus_list[[1]]
for (i in crc_num[-1]) {
  crc_pval_genus <- union(crc_pval_genus,pval_genus_list[[i]])
}
length(crc_pval_genus)

crc_pval_genus_count <- c()
for (i in crc_pval_genus) {
  count <- 0
  for (j in crc_num) {
    if(i %in% pval_genus_list[[j]]) {
      count <- count + 1
    }
  }
  crc_pval_genus_count <- c(crc_pval_genus_count,count)
}

table(crc_pval_genus_count)

count1_crc_pval_genus <- crc_pval_genus[which(crc_pval_genus_count == 1)]
for (i in crc_num) {
  num <- which(crc_num == i)
  pval_samples[num] <- length(pval_genus_list[[i]][pval_genus_list[[i]] %in% count1_crc_pval_genus] )
}


flower_plot <- function(sample, otu_num, core_otu, start, a, b, r, ellipse_col, circle_col) {
  par( bty = 'n', ann = F, xaxt = 'n', yaxt = 'n', mar = c(1,1,1,1))
  plot(c(0,10),c(0,10),type='n')
  n   <- length(sample)
  deg <- 360 / n
  res <- lapply(1:n, function(t){
    draw.ellipse(x = 5 + cos((start + deg * (t - 1)) * pi / 180),
                 y = 5 + sin((start + deg * (t - 1)) * pi / 180),
                 col = ellipse_col[t],
                 border = ellipse_col[t],
                 a = a, b = b, angle = deg * (t - 1))
    text(x = 5 + 2.5 * cos((start + deg * (t - 1)) * pi / 180),
         y = 5 + 2.5 * sin((start + deg * (t - 1)) * pi / 180),
         otu_num[t],cex=2)
    
    if (deg * (t - 1) < 180 && deg * (t - 1) > 0 ) {
      text(x = 5 + 3.3 * cos((start + deg * (t - 1)) * pi / 180),
           y = 5 + 3.3 * sin((start + deg * (t - 1)) * pi / 180),
           sample[t],
           srt = deg * (t - 1) - start,
           adj = 1,
           cex = 2
      )
    } else {
      text(x = 5 + 3.3 * cos((start + deg * (t - 1)) * pi / 180),
           y = 5 + 3.3 * sin((start + deg * (t - 1)) * pi / 180),
           sample[t],
           srt = deg * (t - 1) + start,
           adj = 0,
           cex = 2
      )
    }			
  })
  draw.circle(x = 5, y = 5, r = r, col = circle_col, lwd=2)
}



study.cols <- c("#2ca25f","#8856a7","#e34a33","#3182bd","#F2CC30","#737373","#78c679")

####plot fig1c
flower_plot(sample = paste0("crc",crc_num),
            otu_num = pval_samples,
            core_otu = 43,
            start = 90, a = 1, b = 2,r=1.2,
            ellipse_col = study.cols,circle_col = "white"
)
```




```{r Figure 1e}
zc_samples <- c()
de_samples <- c()
crc_num <- c(1,2,3,4,5,6,7)
for (i in crc_num) {
  crc <- read.table(paste0("d_crc",i,".txt"),header = T,sep="\t")
  zc <- read.table(paste0("h_crc",i,".txt"),header = T,sep="\t")
  meta_crc <- meta_all[which(meta_all$dataset == paste0("crc",i)),]
  de_samples <- c(de_samples,sum(colnames(crc) %in% meta_crc$sample.ID[which(meta_crc$gender == "de")]))
  zc_samples <- c(zc_samples,sum(colnames(zc) %in% meta_crc$sample.ID[which(meta_crc$gender == "zc")]))
  
}
de_samples
zc_samples
##############################################
network_attribute <- function(x) {
  library(igraph)
  net <- graph_from_adjacency_matrix(x,mode="undirected",weighted=T,diag=FALSE)
  igraph.weight <- E(net)$weight
  net <- graph_from_adjacency_matrix(x,mode="undirected",weighted=NULL,diag=FALSE)
  bad.vs = V(net)[degree(net) == 0]
  net = delete.vertices(net, bad.vs)
  #net
  #plot(net,main="Co-occurrence network",vertex.frame.color=NA,vertex.label=NA,edge.width=1,
  #     vertex.size=5,edge.lty=1,edge.curved=TRUE,margin=c(0,0,0,0))
  # network property
  #  The size of the graph (number of edges)
  num.edges = length(E(net)) # length(curve_multiple(igraph))
  #num.edges
  #  Order (number of vertices) of a graph
  num.vertices = length(V(net))# length(diversity(igraph, weights = NULL, vids = V(igraph)))
  #num.vertices
  # (connectance) 
  connectance = edge_density(net,loops=FALSE)# 
  #connectance
  # (Average degree)
  average.degree = mean(igraph::degree(net))# 
  #average.degree
  # Average path length)
  average.path.length = average.path.length(net) # # mean_distance calculates the average path length in a graph
  #average.path.length
  # (Diameter)
  diameter = diameter(net, directed = FALSE, unconnected = TRUE, weights = NULL)
  #diameter
  # edge connectivity / group adhesion
  edge.connectivity = edge_connectivity(net)
  #edge.connectivity
  # (Clustering coefficient)
  clustering.coefficient = transitivity(net) 
  #clustering.coefficient
  no.clusters = no.clusters(net)
  #no.clusters
  # (Betweenness centralization)
  centralization.betweenness = centralization.betweenness(net)$centralization 
  #centralization.betweenness
  # Degree centralization)
  centralization.degree = centralization.degree(net)$centralization
  #centralization.degree
  attribute <- c(num.edges,num.vertices,connectance,average.degree,average.path.length,diameter,edge.connectivity,clustering.coefficient,no.clusters,centralization.betweenness,centralization.degree)
  return(attribute)
}

for (thred in seq(0,0.3,0.03)) {
  case_data_list <- list()
  case_filter_data_list <- list()
  for (m in crc_num) {
    sp_read <- function(x) {
      rownames(x) <- as.character(x$X)
      x <- x[,-1]
      colnames(x) <- rownames(x)
      rownames(x) <- gsub("\\.","-",rownames(x))
      colnames(x) <- gsub("\\.","-",colnames(x))
      rownames(x) <- gsub(" ","-",rownames(x))
      colnames(x) <- gsub(" ","-",colnames(x))
      rownames(x) <- gsub('\\[','',rownames(x))
      colnames(x) <- gsub('\\[','',colnames(x))
      rownames(x) <- gsub('\\]','',rownames(x))
      colnames(x) <- gsub('\\]','',colnames(x))
      return(x)
    }
    sparcc_de <- read.table(paste0("crc",m,"_net.txt"),header = T,sep="\t")
    sp_de <- sp_read(sparcc_de)
    de_filter <- as.matrix(sp_de)
    de_filter[abs(de_filter) < thred] = 0  
    de_filter[de_filter > 0]
    diag(de_filter) <- 1
    de_filter[de_filter > 0] = 1
    de_filter[de_filter < 0] = -1
    case_filter_data_list[[m]] <- de_filter
  }
  crc_network_attribute <- as.data.frame(matrix(nrow=11,ncol=7,NA))
  rownames(crc_network_attribute) <- c("num.edges","num.vertices","connectance","average.degree","average.path.length","diameter","edge.connectivity","clustering.coefficient","no.clusters","centralization.betweenness","centralization.degree")
  colnames(crc_network_attribute) <- paste0("crc",crc_num)
  for (i in crc_num){
    crc_network_attribute[,which(crc_num == i)] = network_attribute(case_filter_data_list[[i]])
  }
  crc_network_attribute[is.na(crc_network_attribute)] = 0
  eval(parse(text = paste("crc_network_attribute_",thred," <- crc_network_attribute",sep="")))
  
}

attribute_names <- c("num.edges","num.vertices","connectance","average.degree","average.path.length","diameter","edge.connectivity","clustering.coefficient","no.clusters","centralization.betweenness","centralization.degree")

attribute_names_vector <- c()
crc_names_vector <- c()
seq_vector <- c()
attribute_num_vector <- c()

for (i in attribute_names){
  for (thred in seq(0,0.3,0.03)) {
    eval(parse(text = paste("crc_network_attribute", "<- crc_network_attribute_",thred,sep="")))
    attribute_num_vector <- c(attribute_num_vector,as.numeric(crc_network_attribute[which(attribute_names == i),]))
    attribute_names_vector <- c(attribute_names_vector,rep(i,7))
    seq_vector <- c(seq_vector,rep(thred,7))
    crc_names_vector <- c(crc_names_vector,names(crc_network_attribute[which(attribute_names == i),]))
  }
}

crc_network_attribute_p <- data.frame(Attributes = attribute_names_vector,Value = attribute_num_vector,Thred = seq_vector,Studies = crc_names_vector)
crc_network_attribute_p$Attributes <- factor(crc_network_attribute_p$Attributes,levels = attribute_names)
crc_network_attribute_p$Studies <- factor(crc_network_attribute_p$Studies)
table(crc_network_attribute_p$Attributes)
crc_network_attribute_p_edge <- crc_network_attribute_p[which(crc_network_attribute_p$Attributes == "num.edges"),]
crc_network_attribute_p_vertex <- crc_network_attribute_p[which(crc_network_attribute_p$Attributes == "num.vertices"),]
crc_network_attribute_p_edge_vertex <- crc_network_attribute_p_edge
names(crc_network_attribute_p_edge_vertex)[2] <- "Edge"
crc_network_attribute_p_edge_vertex$Vertex <- crc_network_attribute_p_vertex$Value
crc_network_attribute_p_degree <- crc_network_attribute_p[which(crc_network_attribute_p$Attributes == "average.degree"),]
crc_network_attribute_p_path <- crc_network_attribute_p[which(crc_network_attribute_p$Attributes == "average.path.length"),]

###plot fig1e
#pdf("fig1f.pdf",width = 6,height = 5)
#edge
ggplot(crc_network_attribute_p_edge_vertex) +
  geom_bar(stat = 'identity',aes(x = Thred, y = Edge, fill = Studies))+
  #geom_line(aes(x=Thred,y=Edge,group=Studies,color=Studies),size=1) +
  #geom_point(aes(x=Thred,y=Edge,group=Studies,size=Vertex,color = Studies),shape=21,stroke =0.5) +
  scale_size_continuous(range = c(5,15)) +
  scale_fill_manual(values =  c("#32955c","#7b4f94","#cc4a33","#3376a8","#e0bf3d","#737373","#72b673"))+
  theme_bw() + 
  theme(
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank(),
    panel.grid.major.y = element_blank(),
    panel.grid.minor.y = element_blank(),)+
    ylab("#Edges") + xlab("Threshold")
#dev.off()
#degree  
#pdf("fig1ff2.pdf",width = 6,height = 5)
ggplot(crc_network_attribute_p_degree) +
  #geom_bar(stat = 'identity',aes(x = Thred, y = Value, fill = Studies))+
  #geom_line(aes(x=Thred,y=Value,group=Studies,color=Studies),size=1) +
  geom_smooth(aes(x=Thred,y=Value,group=Studies,color=Studies),method = 'loess')+
  scale_size_continuous(range = c(5,15)) +
  scale_color_manual(values =  c("#32955c","#7b4f94","#cc4a33","#3376a8","#e0bf3d","#737373","#72b673"))+
  theme_bw() + 
  theme(
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank(),
    panel.grid.major.y = element_blank(),
    panel.grid.minor.y = element_blank(),)+
    ylab("Average degree") + xlab("")
#dev.off()

```
