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

####################################################################
##                        instruction:                            ##
##  put the R markdown files (Figure*.Rmd) and the source data    ##
##  (in the data.zip) in the same directory to run the code       ##
####################################################################

```{r load packages}
library(ggplot2)
library(RColorBrewer)
library(reshape2)
library(ggpubr)
library(MESS)
library(pheatmap)
```



```{r read crc files}
rm(list = ls())
zc_samples <- c()
de_samples <- c()
crc_num <- c(1,2,3,4,5,6,7)
meta_all <- read.csv("ZC.csv",as.is=T) 
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")]))
  
}
case_data_list <- list()
control_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)
  }
  
  #intersect(rownames(x),colnames(x))
  
  sparcc_zc <- read.table(paste0("h_crc",m,"_net.txt"),header = T,sep="\t")
  sparcc_de <- read.table(paste0("d_crc",m,"_net.txt"),header = T,sep="\t")
  
  sp_de <- sp_read(sparcc_de)
  sp_zc <- sp_read(sparcc_zc)
  
  case_data_list[[m]] <- sp_de
  control_data_list[[m]] <- sp_zc
}
crc_num_min <- 1
crc_num_mid <- 4
crc_num_max <- 7
#########################################################
union_genus <- rownames(case_data_list[[1]])
for (i in crc_num) {
  eval(parse(text = paste("union_genus <- union(union_genus,rownames(case_data_list[[",i,"]]))",sep="")))
  
}

union_matrix <- matrix(nrow=length(union_genus),ncol = length(union_genus),0)
rownames(union_matrix) <- union_genus
colnames(union_matrix) <- union_genus
diag(union_matrix) <- 1

case_union_data_list <- list()
for (i in crc_num) {
  fi <- case_data_list[[i]]
  re_union <- union_matrix
  diff_genus <- setdiff(union_genus,rownames(fi))
  rownames(re_union) <- c(rownames(fi),diff_genus)
  colnames(re_union) <- c(rownames(fi),diff_genus)
  for (j in 1:nrow(fi)) {
    re_union[j,1:nrow(fi)] <- as.numeric(fi[j,])
  }
  re_union <- re_union[union_genus,union_genus]
  diag(re_union) <- 1
  case_union_data_list[[i]] <- re_union
  print(i)
}


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

union_matrix <- matrix(nrow=length(union_genus2),ncol = length(union_genus2),0)
rownames(union_matrix) <- union_genus2
colnames(union_matrix) <- union_genus2
diag(union_matrix) <- 1

control_union_data_list <- list()
for (i in crc_num) {
  fi <- control_data_list[[i]]
  re_union <- union_matrix
  diff_genus <- setdiff(union_genus2,rownames(fi))
  rownames(re_union) <- c(rownames(fi),diff_genus)
  colnames(re_union) <- c(rownames(fi),diff_genus)
  for (j in 1:nrow(fi)) {
    re_union[j,1:nrow(fi)] <- as.numeric(fi[j,])
  }
  re_union <- re_union[union_genus2,union_genus2]
  diag(re_union) <- 1
  control_union_data_list[[i]] <- re_union
  print(i)
}

```



```{r Figure 2c}
##case
union_edge_names <- c()
for (i in 1:(nrow(case_union_data_list[[1]])-1)) {
  for (j in (i+1):nrow(case_union_data_list[[1]])) {
    union_edge_names <- c(union_edge_names,paste(union_genus[i],union_genus[j],sep="-"))    
  }
}

union_edge_mat <- matrix(nrow=length(union_edge_names),ncol=7,NA)
rownames(union_edge_mat) <- union_edge_names
colnames(union_edge_mat) <- paste0("crc",crc_num)

for (i in crc_num) {
  union_edge <- c()
  for (k in 1:(nrow(case_union_data_list[[i]])-1)) {
    for (j in (k+1):nrow(case_union_data_list[[i]])) {
      union_edge <- c(union_edge,case_union_data_list[[i]][k,j])    
    }
  }
  union_edge_mat[,which(crc_num ==i)] <- union_edge
}

union_edge_mat <- as.data.frame(union_edge_mat)
dim(union_edge_mat)

union_edge_mat$pool <- NA
for (i in 1:nrow(union_edge_mat)) {
  son <- 0 
  mom <- 0
  for (j in crc_num) {
    v <- (1 - union_edge_mat[i,j]^2)/(de_samples[j] - 1)
    w <- 1/v
    son <- son + w*union_edge_mat[i,j]
    mom <- mom + w
  }
  union_edge_mat$pool[i] <- son/mom
}

##control
union_edge_names <- c()
for (i in 1:(nrow(control_union_data_list[[1]])-1)) {
  for (j in (i+1):nrow(control_union_data_list[[1]])) {
    union_edge_names <- c(union_edge_names,paste(union_genus2[i],union_genus2[j],sep="-"))    
  }
}

zc_edge_mat <- matrix(nrow=length(union_edge_names),ncol=7,NA)
rownames(zc_edge_mat) <- union_edge_names
colnames(zc_edge_mat) <- paste0("zc",crc_num)

for (i in crc_num) {
  union_edge <- c()
  for (k in 1:(nrow(control_union_data_list[[i]])-1)) {
    for (j in (k+1):nrow(control_union_data_list[[i]])) {
      union_edge <- c(union_edge,control_union_data_list[[i]][k,j])    
    }
  }
  zc_edge_mat[,which(crc_num ==i)] <- union_edge
}

zc_edge_mat <- as.data.frame(zc_edge_mat)
dim(zc_edge_mat)

zc_edge_mat$pool <- NA
for (i in 1:nrow(zc_edge_mat)) {
  son <- 0 
  mom <- 0
  for (j in crc_num) {
    v <- (1 - zc_edge_mat[i,j]^2)/(zc_samples[j] - 1)
    w <- 1/v
    son <- son + w*zc_edge_mat[i,j]
    mom <- mom + w
  }
  zc_edge_mat$pool[i] <- son/mom
}

crc.lm <- lm(pool~crc1+crc2+crc3+crc4+crc5+crc6+crc7,data=union_edge_mat)
summary(crc.lm)

zc.lm <- lm(pool~zc1+zc2+zc3+zc4+zc5+zc6+zc7,data=zc_edge_mat)
summary(zc.lm)


crc_error <- as.numeric(summary(crc.lm)$coefficients[,2][-1])
crc_coef <- as.numeric(crc.lm$coefficients[-1])
crc_coef_error_dat <- data.frame(group = rep("crc",7),study = paste("crc",crc_num,sep=""),coef=crc_coef,error=crc_error)

zc_error <- as.numeric(summary(zc.lm)$coefficients[,2][-1])
zc_coef <- as.numeric(zc.lm$coefficients[-1])
zc_coef_error_dat <- data.frame(group = rep("zc",7),study = paste("crc",crc_num,sep=""),coef=zc_coef,error=zc_error)

coef_error_dat <- rbind(crc_coef_error_dat,zc_coef_error_dat)
coef_error_dat$study <- as.factor(coef_error_dat$study)

##plot figure 2c
#pdf("fig2c.pdf")
ggplot(coef_error_dat,aes(x=study,y=coef,group=group,fill=group)) +
  geom_bar(stat = 'identity')+
  facet_grid(group~.,scales = "free")+
  scale_fill_manual(values =  c("#d9544d","#C2C2C2"), label = c("case","control"))+
  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(),
    axis.title.x=element_blank(),
    axis.title.y = element_blank())
#dev.off()
```
```{r Figure 2d}
vi = read.table("univariate_dis_case.txt",header = T,sep = '\t',row.names = 1)
sp = read.table("simple_dis_case.txt",header = T,sep = '\t',row.names = 1)
vi.nmds = read.table("univariate_nmds.txt",header = T,sep = '\t',row.names = 1)
sp.nmds = read.table("simple_nmds.txt",header = T,sep = '\t',row.names = 1)
colnames(vi) = rownames(vi)
colnames(sp) = rownames(sp)
aa = data.frame(vi.nmds[,4])
aa$vi.nmds...4. = factor(aa$vi.nmds...4.)

##plot fig2d
#univariate
pheatmap(vi,cluster_rows = F,cluster_cols = F,
         annotation_row = aa,annotation_col = aa,
         show_rownames = F,show_colnames = F)
#unweighted
pheatmap(sp,cluster_rows = F,cluster_cols = F,
         annotation_row = aa,annotation_col = aa,
         show_rownames = F,show_colnames = F)
```



```{r Figure 2e}
vi = read.table("univariate_dis_case.txt",header = T,sep = '\t',row.names = 1)
sp = read.table("simple_dis_case.txt",header = T,sep = '\t',row.names = 1)
vi.nmds = read.table("univariate_nmds.txt",header = T,sep = '\t',row.names = 1)
sp.nmds = read.table("simple_nmds.txt",header = T,sep = '\t',row.names = 1)
colnames(vi) = rownames(vi)
colnames(sp) = rownames(sp)
bb = vi
bb$name = rownames(bb)
dd = melt(bb,id.vars = "name")
dd$type1 = vi.nmds[as.character(dd$name),4]
dd$type2 = vi.nmds[as.character(dd$variable),4]
vi.dis = dd[which(dd$type1 == dd$type2),]
bb = sp
bb$name = rownames(bb)
dd = melt(bb,id.vars = "name")
dd$type1 = vi.nmds[as.character(dd$name),4]
dd$type2 = vi.nmds[as.character(dd$variable),4]
sp.dis = dd[which(dd$type1 == dd$type2),]
vi.dis$class = "univariate weighting"
sp.dis$class = "unweighted"
mydata = rbind(vi.dis,sp.dis)

##plot fig2e
ggplot(mydata[-which(mydata$type2 == 7),])+
  geom_boxplot(aes(factor(type2),value,fill = class))+
  theme_bw()+
  scale_fill_manual(values = c("#C3514A","#C2C1C0"))
  theme(panel.grid = element_blank())

```


```{r Figure 2f}
dis.uni = read.table("cor_univariate.txt",header = T,sep = '\t')
dis.sim = read.table("cor_simple.txt",header = T,sep = '\t')
dis.uni$type = "univariate weighting"
dis.sim$type = "unweighted"
dis.all = rbind(dis.uni,dis.sim)
ggplot(dis.all,aes(sampleDis,distance,color = type))+
  geom_smooth(method = 'loess')+
  scale_color_manual(values = c("#C3514A","#C2C1C0"))+
  theme_bw()+theme(panel.grid = element_blank())+
  labs(x = "Sample dissimilarity", y = "Network distance")
#cor.test(dis.uni$distance,dis.uni$sampleDis)
#cor.test(dis.sim$distance,dis.sim$sampleDis)
```


```{r Figure 2g}
net = read.table(paste0("d_netmoss_nor_result.txt"),header = T,sep = '\t')
combat = read.table(paste0("d_combat_nor_result.txt"),header = T,sep = '\t')
limma = read.table(paste0("d_limma_nor_result.txt"),header = T,sep = '\t')
merge = read.table(paste0("d_merge_nor_result.txt"),header = T,sep = '\t')

##crc
case_name = "crc"
crc_num <- c(1,2,3,4,5,6,7)
control_samples <- c()
case_samples <- c()
for (i in crc_num) {
  crc <- read.table(paste0("d_",case_name,i,".txt"),header = T,sep="\t")
  control <- read.table(paste0("h_",case_name,i,".txt"),header = T,sep="\t")
  case_samples <- c(case_samples,length(colnames(crc))-1)
  control_samples <- c(control_samples,length(colnames(control))-1)
  
}
case_data_list <- list()
control_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_control <- read.table(paste0("h_",case_name,m,"_net2.txt"),header = T,sep="\t")
  sparcc_case <- read.table(paste0("d_",case_name,m,"_net2.txt"),header = T,sep="\t")
  sp_case <- sp_read(sparcc_case)
  sp_control <- sp_read(sparcc_control)
  case_data_list[[m]] <- sp_case
  control_data_list[[m]] <- sp_control
}

net = sp_read(net)
combat = sp_read(combat)
merge = sp_read(merge)
limma = sp_read(limma)


union_genus <- rownames(case_data_list[[1]])
for (i in crc_num) {
  eval(parse(text = paste("union_genus <- union(union_genus,rownames(case_data_list[[",i,"]]))",sep="")))
  
}
union_matrix <- matrix(nrow=length(union_genus),ncol = length(union_genus),0)
rownames(union_matrix) <- union_genus
colnames(union_matrix) <- union_genus
diag(union_matrix) <- 1

case_union_data_list <- list()
for (i in crc_num) {
  fi <- case_data_list[[i]]
  re_union <- union_matrix
  diff_genus <- setdiff(union_genus,rownames(fi))
  rownames(re_union) <- c(rownames(fi),diff_genus)
  colnames(re_union) <- c(rownames(fi),diff_genus)
  for (j in 1:nrow(fi)) {
    re_union[j,1:nrow(fi)] <- as.numeric(fi[j,])
  }
  re_union <- re_union[union_genus,union_genus]
  diag(re_union) <- 1
  case_union_data_list[[i]] <- re_union
  print(i)
}

union_genus2 <- rownames(control_data_list[[1]])
for (i in crc_num) {
  eval(parse(text = paste("union_genus2 <- union(union_genus2,rownames(control_data_list[[",i,"]]))",sep="")))
  
}
union_matrix <- matrix(nrow=length(union_genus2),ncol = length(union_genus2),0)
rownames(union_matrix) <- union_genus2
colnames(union_matrix) <- union_genus2
diag(union_matrix) <- 1

control_union_data_list <- list()
for (i in crc_num) {
  fi <- control_data_list[[i]]
  re_union <- union_matrix
  diff_genus <- setdiff(union_genus2,rownames(fi))
  rownames(re_union) <- c(rownames(fi),diff_genus)
  colnames(re_union) <- c(rownames(fi),diff_genus)
  for (j in 1:nrow(fi)) {
    re_union[j,1:nrow(fi)] <- as.numeric(fi[j,])
  }
  re_union <- re_union[union_genus2,union_genus2]
  diag(re_union) <- 1
  control_union_data_list[[i]] <- re_union
  print(i)
}

genus.all1 = intersect(rownames(net),rownames(combat))
genus.all2 = intersect(rownames(merge),rownames(limma))
genus.all3 = intersect(genus.all1,genus.all2)
genus.all4 = intersect(union_genus,union_genus2)
genus.all = intersect(genus.all3,genus.all4)

##distance
case_union_data_list[[1]] = case_union_data_list[[1]][as.character(rownames(case_union_data_list[[1]])) %in% as.character(genus.all),
                                                      as.character(colnames(case_union_data_list[[1]])) %in% as.character(genus.all)]
case_union_data_list[[2]] = case_union_data_list[[2]][as.character(rownames(case_union_data_list[[2]])) %in% as.character(genus.all),
                                                      as.character(colnames(case_union_data_list[[2]])) %in% as.character(genus.all)]
case_union_data_list[[3]] = case_union_data_list[[3]][as.character(rownames(case_union_data_list[[3]])) %in% as.character(genus.all),
                                                      as.character(colnames(case_union_data_list[[3]])) %in% as.character(genus.all)]
case_union_data_list[[4]] = case_union_data_list[[4]][as.character(rownames(case_union_data_list[[4]])) %in% as.character(genus.all),
                                                      as.character(colnames(case_union_data_list[[4]])) %in% as.character(genus.all)]
case_union_data_list[[5]] = case_union_data_list[[5]][as.character(rownames(case_union_data_list[[5]])) %in% as.character(genus.all),
                                                      as.character(colnames(case_union_data_list[[5]])) %in% as.character(genus.all)]
case_union_data_list[[6]] = case_union_data_list[[6]][as.character(rownames(case_union_data_list[[6]])) %in% as.character(genus.all),
                                                      as.character(colnames(case_union_data_list[[6]])) %in% as.character(genus.all)]
case_union_data_list[[7]] = case_union_data_list[[7]][as.character(rownames(case_union_data_list[[7]])) %in% as.character(genus.all),
                                                      as.character(colnames(case_union_data_list[[7]])) %in% as.character(genus.all)]
control_union_data_list[[1]] = control_union_data_list[[1]][as.character(rownames(control_union_data_list[[1]])) %in% as.character(genus.all),
                                                            as.character(colnames(control_union_data_list[[1]])) %in% as.character(genus.all)]
control_union_data_list[[2]] = control_union_data_list[[2]][as.character(rownames(control_union_data_list[[2]])) %in% as.character(genus.all),
                                                            as.character(colnames(control_union_data_list[[2]])) %in% as.character(genus.all)]
control_union_data_list[[3]] = control_union_data_list[[3]][as.character(rownames(control_union_data_list[[3]])) %in% as.character(genus.all),
                                                            as.character(colnames(control_union_data_list[[3]])) %in% as.character(genus.all)]
control_union_data_list[[4]] = control_union_data_list[[4]][as.character(rownames(control_union_data_list[[4]])) %in% as.character(genus.all),
                                                            as.character(colnames(control_union_data_list[[4]])) %in% as.character(genus.all)]
control_union_data_list[[5]] = control_union_data_list[[5]][as.character(rownames(control_union_data_list[[5]])) %in% as.character(genus.all),
                                                            as.character(colnames(control_union_data_list[[5]])) %in% as.character(genus.all)]
control_union_data_list[[6]] = control_union_data_list[[6]][as.character(rownames(control_union_data_list[[6]])) %in% as.character(genus.all),
                                                            as.character(colnames(control_union_data_list[[6]])) %in% as.character(genus.all)]
control_union_data_list[[7]] = control_union_data_list[[7]][as.character(rownames(control_union_data_list[[7]])) %in% as.character(genus.all),
                                                            as.character(colnames(control_union_data_list[[7]])) %in% as.character(genus.all)]
net = net[as.character(rownames(net)) %in% as.character(genus.all),as.character(colnames(net)) %in% as.character(genus.all)]
combat = combat[as.character(rownames(combat)) %in% as.character(genus.all),as.character(colnames(combat)) %in% as.character(genus.all)]
limma = limma[as.character(rownames(limma)) %in% as.character(genus.all),as.character(colnames(limma)) %in% as.character(genus.all)]
merge = merge[as.character(rownames(merge)) %in% as.character(genus.all),as.character(colnames(merge)) %in% as.character(genus.all)]

data.dis = data.frame()
for (k in 1:7)
{
  data.dis[1,k] = cmd(as.matrix(case_union_data_list[[k]]),as.matrix(net))
  data.dis[2,k] = cmd(as.matrix(case_union_data_list[[k]]),as.matrix(combat))
  data.dis[3,k] = cmd(as.matrix(case_union_data_list[[k]]),as.matrix(limma))
  data.dis[4,k] = cmd(as.matrix(case_union_data_list[[k]]),as.matrix(merge))
}
colnames(data.dis) = c("crc1","crc2","crc3","crc4","crc5","crc6","crc7")
rownames(data.dis) = c("net.cor","combat.cor","limma.cor","merge.cor")
data.dis$name = rownames(data.dis)
data.dis2 = melt(data.dis,id.vars = "name")
data.dis2$value2 = 1 - data.dis2$value
data.dis2$value2 = (data.dis2$value2-min(data.dis2$value2))/(max(data.dis2$value2)-min(data.dis2$value2))
data.dis2$x = as.numeric(substr(data.dis2$variable,4,5))


###plot fig2g#
ggplot(data.dis2,aes(factor(x),value,color = name,group = name))+
  geom_point(position = position_jitter(.2))+
  geom_line()+
  theme_bw()+
  theme(panel.grid = element_blank())+
  scale_color_brewer(palette = "Set1")
```