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

```{r load packages}
library(reshape2)
library(ggplot2)
library(MESS)
library(plotROC)
library(WeightedROC)
```



```{r read markers}
rm(list = ls())
##set threshold
threshold = 0.12
#maker
gutmdisorder_crc <- read.csv("data.csv",header = F,stringsAsFactors = F)
table(gutmdisorder_crc$V9[which(gutmdisorder_crc$V6 == "genus")])
target_genus <- gutmdisorder_crc$V5[which(gutmdisorder_crc$V6 == "genus")]
target_g <- c()
for (i in 1:length(target_genus)) {
  target_g <- c(target_g,strsplit(target_genus," ")[[i]][1])
}
target_g <- unique(target_g)
target_gg <- paste("g__",target_g,sep="")
```


```{r Figure 4b}
m.case = read.table("case_class.txt",header = T,sep = '\t')
m.control = read.table("control_class.txt",header = T,sep = '\t')
m.case1 = dcast(m.case[,c(4,3,2)],mod~phylum,mean,fill = 0)
m.case2 = melt(m.case1,id.vars = "mod")
##plot figure 4b case
ggplot(m.case2,aes(mod,value,fill = variable))+
  geom_bar(stat = 'identity',position = 'fill',width = 1,color = "black")+
  coord_polar(theta = "y")+
  theme_bw()+
  theme(panel.grid = element_blank())+
  theme(panel.border = element_blank())+
  theme(axis.text = element_blank(),axis.ticks = element_blank())+
  scale_fill_manual(values=c("#ffffb3","#bebada","#8dd3c7","#fb8072","#80b1d3"))+
  labs(x = "", y = "case")
m.control1 = dcast(m.control[,c(4,3,2)],mod~phylum,mean,fill = 0)
m.control2 = melt(m.control1,id.vars = "mod")
##plot figure 4b control
ggplot(m.control2,aes(mod,value,fill = variable))+
  geom_bar(stat = 'identity',position = 'fill',width = 1,color = "black")+
  coord_polar(theta = "y")+
  theme_bw()+
  theme(panel.grid = element_blank())+
  theme(panel.border = element_blank())+
  scale_fill_manual(values=c("#ffffb3","#bebada","#8dd3c7","#fb8072","#80b1d3")) +
  theme(axis.text = element_blank(),axis.ticks = element_blank())+
  labs(x = "", y = "control")
```



```{r Figure 4c}
nodes_result = read.table("node_result_netmoss.txt",header = T,sep = '\t')
nodes_result$target = 0
rownames(nodes_result) = nodes_result$taxon_names
nodes_result[as.character(rownames(nodes_result)) %in% as.character(target_gg),"target"] = 1
##degree
case.sample = read.table("d_netmoss_nor_result.txt",header = T,sep = '\t',row.names = 1)
case.node = data.frame()
for (rr in 1:nrow(case.sample))
{
  case.node[rr,1] = rownames(case.sample)[rr]
  case.node[rr,2] = length(which(case.sample[rr,] != 0))
}
control.sample = read.table("h_netmoss_nor_result.txt",header = T,sep = '\t',row.names = 1)
control.node = data.frame()
for (rr in 1:nrow(control.sample))
{
  control.node[rr,1] = rownames(control.sample)[rr]
  control.node[rr,2] = length(which(control.sample[rr,] != 0))
}
case.node$type = "case"
colnames(case.node)[1] = "node"
colnames(case.node)[2] = "degree"
rownames(case.node) = case.node$node
case.node[as.character(target_gg),"target"] = 1
case.node[which(is.na(case.node$target)),"target"] = 0
case.node = case.node[-which(is.na(case.node[,1])),]
control.node$type = "control"
colnames(control.node)[1] = "node"
colnames(control.node)[2] = "degree"
rownames(control.node) = control.node$node
control.node[as.character(target_gg),"target"] = 1
control.node[which(is.na(control.node$target)),"target"] = 0
control.node = control.node[-which(is.na(control.node[,1])),]
crc_node = rbind(case.node,control.node)
##plot figure 4c
ggplot(crc_node,aes(type, degree, fill = factor(target)))+
  geom_boxplot()+
  theme_bw()+
  theme(panel.grid = element_blank())+
  labs(x = "", y = "Degree")+
  scale_fill_manual(values=c("#C2C2C2","#d9544d"),label = c("nonmarker","marker"))
#test
#case1 = crc_node[which(crc_node$type == "case"),]
#control1 = crc_node[which(crc_node$type == "control"),]
#wilcox.test(as.numeric(as.character(case1[which(case1$target == 0),"degree"])),
#            as.numeric(as.character(case1[which(case1$target == 1),"degree"])))
#wilcox.test(as.numeric(as.character(control1[which(control1$target == 0),"degree"])),
#            as.numeric(as.character(control1[which(control1$target == 1),"degree"])))


```


```{r Figure 4d}
crc = read.table("node_result_netmoss.txt",header = T,sep = '\t',row.names = 1)
crc2 = crc[,c(11,12,13)]
crc2$target = 0
crc2$net = 0
crc2$ab = 0
crc2[as.character(rownames(crc)) %in% as.character(target_gg),"target"] = 1
crc2[which(crc2$netzgo.scale > threshold),"net"] = 2
crc2[which(crc2$padj < 0.05),"ab"] = 4
crc2$sum = crc2$target+crc2$net+crc2$ab
crc2 = crc2[order(crc2$netzgo.scale,decreasing = T),]

##plot figure 4d density
#NetMoss
ggplot(crc2,aes(x = netzgo.scale))+geom_density(fill = "#B0C9DE")+
  theme_bw()+
  theme(panel.grid = element_blank())+
  labs(x = "NetMoss score", y = "Density")
#Wilcoxon
ggplot(crc2,aes(x = sqrt.log.padj))+geom_density(fill = "#C9E3C1")+
  theme_bw()+
  labs(x = "Sqrt(-log10(q))", y = "Density")+
  theme(panel.grid = element_blank())

##plot figure 4d point
#pdf("fig4d_point.pdf",width = 6,height = 5)
ggplot(crc2,aes(x=netzgo.scale,y=sqrt.log.padj,color=factor(target)))+
  geom_point(size = 2.5) +
  theme_bw() +
  theme(axis.text.x=element_text(colour="black",size=14), 
        axis.text.y=element_text(size=14), 
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  geom_vline(aes(xintercept= threshold), colour="#0663A4", linetype="dashed") +
  geom_hline(aes(yintercept= sqrt(-log10(0.05))), colour="blue", linetype="dashed")+
  scale_y_continuous(expand = c(0,0))+
  scale_x_continuous(expand = c(0,0))+
  scale_color_manual(values = c("#c2c1c1","#cd524b"), label = c("nonmarker","marker"))+
  labs(x = "NetMoss score", y = "Sqrt(-log10(q))")
#dev.off()

target = crc2[which(crc2$target == 1),]
y2 = data.frame()
x = 0.01
for (k in 1:100)
{
  y2[k,1] = x
  y2[k,2] = nrow(target[which(target$netzgo.scale > x),])
  y2[k,3] = nrow(crc2[which(crc2$netzgo.scale > x),])-nrow(target[which(target$netzgo.scale > x),])
  y2[k,4] = nrow(target)-nrow(target[which(target$netzgo.scale > x),])
  y2[k,5] = nrow(crc2)-nrow(crc2[which(crc2$netzgo.scale > x),])-nrow(target)+nrow(target[which(target$netzgo.scale > x),])
  y2[k,6] = y2[k,2]/(y2[k,2]+y2[k,4])
  y2[k,7] = y2[k,5]/(y2[k,3]+y2[k,5])
  x = x + 0.01
}
colnames(y2) = c("threshold","TP","FP","FN","TN","sensitivity","specificity")
y3 = melt(y2[,c(1,6,7)],id.vars = "threshold")

##plot fig4d line#################### 
ggplot(y3,aes(threshold,value,group = variable, color = variable))+
  theme_bw()+
  theme(panel.grid = element_blank())+
  labs(x = "NetMoss",y = "Value")+
  geom_smooth(method = 'loess')+
  scale_color_manual(values = c("#DE9528","#40B389"))
```


```{r Figure 4e}
crc2$class = "4None"
crc2[which(crc2$sum == 4),"class"] = "3Statistical"
crc2[which(crc2$sum == 5),"class"] = "3Statistical"
crc2[which(crc2$sum == 2),"class"] = "1NetmMoss"
crc2[which(crc2$sum == 3),"class"] = "1NetmMoss"
crc2[which(crc2$sum == 6),"class"] = "2Both"
crc2[which(crc2$sum == 7),"class"] = "2Both"
crc2$genus = rownames(crc2)
crc2$genus = factor(crc2$genus,levels = rev(crc2$genus))

crc3 = crc2[which(crc2$target == 1),]
crc3 = crc3[order(crc3$class),]
crc3$genus = factor(crc3$genus,levels = rev(crc3$genus))
crc3 = crc3[-which(crc3$class == "4None"),]

##plot figure 4e
#pdf("fig4e.pdf",width = 6,height = 8)
ggplot(crc3,aes(genus,netzgo.scale,color = class))+
  geom_point(shape = 18,size = 2)+
  theme_bw()+
  theme(panel.grid = element_blank())+
  coord_flip()+
  geom_hline(yintercept = threshold)
#dev.off()
```


```{r Figure 4f}
crc1 = read.table("netmoss_roc.txt",header = T,sep = '\t')
crc2 = read.table("rf_roc.txt",header = T,sep = '\t')
crc3 = read.table("plsda_roc.txt",header = T,sep = '\t')
crc4 = read.table("merge_roc.txt",header = T,sep = '\t')
crc5 = read.table("limma_roc.txt",header = T,sep = '\t')
crc6 = read.table("combat_roc.txt",header = T,sep = '\t')
crc1$type = "NetMoss"
crc2$type = "RandomForest"
crc3$type = "PLS-DA"
crc4$type = "Wilcoxon"
crc5$type = "limma"
crc6$type = "combat"
mydata = rbind(crc1,crc2,crc3,crc4,crc5,crc6)

###auc
auc.all = data.frame(type = c("NetMoss","Random Forest","PLS-DA","WilCoxon","limma","combat"),
                     auc = c(WeightedAUC(crc1),
                             WeightedAUC(crc2),
                             WeightedAUC(crc3),
                             WeightedAUC(crc4),
                             WeightedAUC(crc5),
                             WeightedAUC(crc6)))
###plot figure 4f
ggplot()+ 
  geom_path(aes(FPR, TPR,color = type), data=mydata)+
  coord_equal()+
  style_roc()+
  annotate("text",x = .75, y = .30, ## 
           label = paste("AUC of NetMoss = ", round(auc.all[1,2], 2)))+
  annotate("text",x = .75, y = .25, ##
           label = paste("AUC of Random Forest = ", round(auc.all[2,2], 2)))+
  annotate("text",x = .75, y = .20, ## 
           label = paste("AUC of PLS-DA = ", round(auc.all[3,2], 2)))+
  annotate("text",x = .75, y = .15, ## 
           label = paste("AUC of Wilcoxon = ", round(auc.all[4,2], 2)))+
  annotate("text",x = .75, y = .10, ## 
           label = paste("AUC of limma = ", round(auc.all[5,2], 2)))+
  annotate("text",x = .75, y = .05, ## 
           label = paste("AUC of combat = ", round(auc.all[6,2], 2)))
```


```{r Figure 4g-4i}
##netmoss
crc1 = read.table("netmoss_roc_crc1.txt",header = T,sep = '\t')
crc2 = read.table("netmoss_roc_crc2.txt",header = T,sep = '\t')
crc3 = read.table("netmoss_roc_crc3.txt",header = T,sep = '\t')
crc4 = read.table("netmoss_roc_crc4.txt",header = T,sep = '\t')
crc5 = read.table("netmoss_roc_crc5.txt",header = T,sep = '\t')
crc6 = read.table("netmoss_roc_crc6.txt",header = T,sep = '\t')
crc7 = read.table("netmoss_roc_crc7.txt",header = T,sep = '\t')

crc1$type = "CRC1"
crc2$type = "CRC2"
crc3$type = "CRC3"
crc4$type = "CRC4"
crc5$type = "CRC5"
crc6$type = "CRC6"
crc7$type = "CRC7"
mydata = rbind(crc1,crc2,crc3,crc4,crc5,crc6,crc7)

##plot figure 4g
p = ggplot(mydata, aes(d = test_orig.select, m = prob, color = type)) + 
  geom_roc(n.cuts = 0) + 
  style_roc()
auc<-calc_auc(p)
p+annotate("text",x = .75, y = .30, ## 注释text的位置
           label = paste("AUC of CRC1 = ", round(calc_auc(p)$AUC[1], 2)))+
  annotate("text",x = .75, y = .25, ## 注释text的位置
           label = paste("AUC of CRC2 = ", round(calc_auc(p)$AUC[2], 2)))+
  annotate("text",x = .75, y = .20, ## 注释text的位置
           label = paste("AUC of CRC3 = ", round(calc_auc(p)$AUC[3], 2)))+
  annotate("text",x = .75, y = .15, ## 注释text的位置
           label = paste("AUC of CRC4 = ", round(calc_auc(p)$AUC[4], 2)))+
  annotate("text",x = .75, y = .10, ## 注释text的位置
           label = paste("AUC of CRC5 = ", round(calc_auc(p)$AUC[5], 2)))+
  annotate("text",x = .75, y = .05, ## 注释text的位置
           label = paste("AUC of CRC6 = ", round(calc_auc(p)$AUC[6], 2)))+
  annotate("text",x = .75, y = 0, ## 注释text的位置
           label = paste("AUC of CRC7 = ", round(calc_auc(p)$AUC[7], 2)))

##plsda
crc1 = read.table("plsda_roc_crc1.txt",header = T,sep = '\t')
crc2 = read.table("plsda_roc_crc2.txt",header = T,sep = '\t')
crc3 = read.table("plsda_roc_crc3.txt",header = T,sep = '\t')
crc4 = read.table("plsda_roc_crc4.txt",header = T,sep = '\t')
crc5 = read.table("plsda_roc_crc5.txt",header = T,sep = '\t')
crc6 = read.table("plsda_roc_crc6.txt",header = T,sep = '\t')
crc7 = read.table("plsda_roc_crc7.txt",header = T,sep = '\t')

crc1$type = "CRC1"
crc2$type = "CRC2"
crc3$type = "CRC3"
crc4$type = "CRC4"
crc5$type = "CRC5"
crc6$type = "CRC6"
crc7$type = "CRC7"
mydata = rbind(crc1,crc2,crc3,crc4,crc5,crc6,crc7)

##plot figure 4h
p = ggplot(mydata, aes(d = test_orig.select, m = prob, color = type)) + 
  geom_roc(n.cuts = 0) + 
  style_roc()
auc<-calc_auc(p)
p+annotate("text",x = .75, y = .30, ## 注释text的位置
           label = paste("AUC of CRC1 = ", round(calc_auc(p)$AUC[1], 2)))+
  annotate("text",x = .75, y = .25, ## 注释text的位置
           label = paste("AUC of CRC2 = ", round(calc_auc(p)$AUC[2], 2)))+
  annotate("text",x = .75, y = .20, ## 注释text的位置
           label = paste("AUC of CRC3 = ", round(calc_auc(p)$AUC[3], 2)))+
  annotate("text",x = .75, y = .15, ## 注释text的位置
           label = paste("AUC of CRC4 = ", round(calc_auc(p)$AUC[4], 2)))+
  annotate("text",x = .75, y = .10, ## 注释text的位置
           label = paste("AUC of CRC5 = ", round(calc_auc(p)$AUC[5], 2)))+
  annotate("text",x = .75, y = .05, ## 注释text的位置
           label = paste("AUC of CRC6 = ", round(calc_auc(p)$AUC[6], 2)))+
  annotate("text",x = .75, y = 0, ## 注释text的位置
           label = paste("AUC of CRC7 = ", round(calc_auc(p)$AUC[7], 2)))


##random forest
crc1 = read.table("rf_roc_crc1.txt",header = T,sep = '\t')
crc2 = read.table("rf_roc_crc2.txt",header = T,sep = '\t')
crc3 = read.table("rf_roc_crc3.txt",header = T,sep = '\t')
crc4 = read.table("rf_roc_crc4.txt",header = T,sep = '\t')
crc5 = read.table("rf_roc_crc5.txt",header = T,sep = '\t')
crc6 = read.table("rf_roc_crc6.txt",header = T,sep = '\t')
crc7 = read.table("rf_roc_crc7.txt",header = T,sep = '\t')

crc1$type = "CRC1"
crc2$type = "CRC2"
crc3$type = "CRC3"
crc4$type = "CRC4"
crc5$type = "CRC5"
crc6$type = "CRC6"
crc7$type = "CRC7"
mydata = rbind(crc1,crc2,crc3,crc4,crc5,crc6,crc7)

##plot figure 4i
p = ggplot(mydata, aes(d = test_orig.select, m = prob, color = type)) + 
  geom_roc(n.cuts = 0) + 
  style_roc()
auc<-calc_auc(p)
p+annotate("text",x = .75, y = .30, ## 注释text的位置
           label = paste("AUC of CRC1 = ", round(calc_auc(p)$AUC[1], 2)))+
  annotate("text",x = .75, y = .25, ## 注释text的位置
           label = paste("AUC of CRC2 = ", round(calc_auc(p)$AUC[2], 2)))+
  annotate("text",x = .75, y = .20, ## 注释text的位置
           label = paste("AUC of CRC3 = ", round(calc_auc(p)$AUC[3], 2)))+
  annotate("text",x = .75, y = .15, ## 注释text的位置
           label = paste("AUC of CRC4 = ", round(calc_auc(p)$AUC[4], 2)))+
  annotate("text",x = .75, y = .10, ## 注释text的位置
           label = paste("AUC of CRC5 = ", round(calc_auc(p)$AUC[5], 2)))+
  annotate("text",x = .75, y = .05, ## 注释text的位置
           label = paste("AUC of CRC6 = ", round(calc_auc(p)$AUC[6], 2)))+
  annotate("text",x = .75, y = 0, ## 注释text的位置
           label = paste("AUC of CRC7 = ", round(calc_auc(p)$AUC[7], 2)))

```

