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

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

```



```{r FIgure 5a}
rm(list = ls())
a = read.table("disease_all.txt",header = T,sep = '\t')
tax = read.table("classification.txt",header = T,sep = '\t')
b = a[,c(1,12,13,15)]
b[,2] = as.numeric(as.character(b[,2]))
b[,3] = as.numeric(as.character(b[,3]))
b$class1 = 0
b$class2 = 0
b[which(b$netzgo.scale > 0.3),"class1"] = 1
b[which(b$padj < 0.01),"class2"] = 2
b$class3 = b$class2+b$class1
dd = b[-which(b$class3 == 0),]
dd[which(dd$class3 == 1),"class"] = "diff-net"
dd[which(dd$class3 == 2),"class"] = "diff-abd"
dd[which(dd$class3 == 3),"class"] = "diff-both"
data.plot = dcast(dd[,c(1,4,7)],type~taxon_names,mean,fill = 0)
rownames(data.plot) = data.plot[,1]
data.plot = data.plot[,-1]

##ordered by type
data1 = data.frame(t(data.plot))
data1$sum = rowSums(data1)
data1 = data1[order(data1$sum,decreasing = T),]
data2 = data.frame(t(data1[,-ncol(data1)]))

##ordered by family
anno = tax[as.character(tax$genus) %in% as.character(colnames(data2)),]
anno = anno[!duplicated(anno$genus),]
anno_col = data.frame(anno[,5])
colnames(anno_col)[1] = "family"
rownames(anno_col) = anno$genus

anno_col[which(anno_col$family == "f__Lachnospiraceae"),"type"] = "Lachnospiraceae"
anno_col[which(anno_col$family == "f__Ruminococcaceae"),"type"] = "Ruminococcaceae"
anno_col[which(anno_col$family == "f__Enterobacteriaceae"),"type"] = "Enterobacteriaceae"
anno_col[which(is.na(anno_col$type)),"type"] = "Others"

anno_col = data.frame(anno_col[,2])
colnames(anno_col)[1] = "family"
rownames(anno_col) = anno$genus
data3 = data.frame(t(data2))
data3$family = anno_col[as.character(rownames(data3)),1]
data3[which(is.na(data3$family)),"family"] = "Others"

data4 = data.frame(t(data3[,-ncol(data3)]))
my.color = colorRampPalette(brewer.pal( 8, "Spectral"))(length(unique(anno_col$family)))
names(my.color) = unique(data3$family)
ann_color = list(c(my.color))
names(ann_color) = "family"

##plot figure5a
pheatmap(data2,color = c("white","#B0C9DE","#C9E3C1","#F4B1AB"),
         cluster_rows = T,cluster_cols = F,
         show_colnames = F,annotation_col = anno_col,
         annotation_colors = ann_color)
```



```{r Figure 5b}
mydata1 = read.table("5disease_ratio.txt",header = T,sep = '\t')

##plot figure 5b
ggplot(mydata1,aes(disease,ratio,fill = disease))+
  geom_bar(stat = 'identity')+
  labs(x = "", y = "Multi-related bacteria rate")+
  theme_bw()+
  theme(panel.grid = element_blank())+
  guides(fill = F)+
  scale_fill_manual(values = c("#984EA3","#E41A1C","#377EB8","#4DAF4A","#FF7F00"))
```



```{r Figure 5d}
case = read.table("node_case.txt",header = T,sep = '\t')
control = read.table("node_control.txt",header = T,sep = '\t')
case$state = "case"
control$state = "control"
mydata = rbind(case,control)

#degree box
my_comparisons = list(c("specific","multiple"))

##plot figure 5d
ggplot(mydata[which(mydata$degree != 0),],aes(group,degree,fill = group))+
  geom_boxplot()+
  facet_grid(.~state)+
  stat_compare_means(comparisons = my_comparisons)+
  theme_bw()+
  theme(panel.grid = element_blank())+
  labs(x = "", y = "Degree")+
  scale_fill_manual(values = c("#CD524B","#C2C1C1"))
```


```{r Figure 5e}
aids = read.table("node_result_HIV.txt",header = T,sep = '\t')
asd = read.table("node_result_ASD.txt",header = T,sep = '\t')
crc = read.table("node_result_CRC.txt",header = T,sep = '\t')
cr = read.table("node_result_CR.txt",header = T,sep = '\t')
ob = read.table("node_result_OB.txt",header = T,sep = '\t')
case = read.table("node_case.txt",header = T,sep = '\t')
control = read.table("node_control.txt",header = T,sep = '\t')


###multiple
case.m = case[which(case$group == "multiple"),]
control.m = control[which(control$group == "multiple"),]
case.m2 = dcast(case.m[,c(1,2,4)],genus~type,sum,fill = 0)
control.m2 = dcast(control.m[,c(1,2,4)],genus~type,sum,fill = 0)
gg = intersect(as.character(case.m2$genus),as.character(control.m2$genus))

ob1 = ob[as.character(ob$taxon_names) %in% as.character(gg),c(1,12)]
cr1 = cr[as.character(cr$taxon_names) %in% as.character(gg),c(1,12)]
crc1 = crc[as.character(crc$taxon_names) %in% as.character(gg),c(1,12)]
asd1 = asd[as.character(asd$taxon_names) %in% as.character(gg),c(1,12)]
aids1 = aids[as.character(aids$taxon_names) %in% as.character(gg),c(1,12)]

ob1$type = "OB"
cr1$type = "CR"
crc1$type = "CRC"
asd1$type = "ASD"
aids1$type = "AIDS"

mydata = rbind(ob1,asd1,aids1,cr1,crc1)
mydata[which(is.na(mydata$netzgo.scale)),"netzgo.scale"] = 0
mydata2 = dcast(mydata[,c(1,3,2)],taxon_names~type,sum,fill = 0)
mydata.multi = mydata2


#######single
case.m = case[which(case$group != "multiple"),]
control.m = control[which(control$group != "multiple"),]
case.m2 = dcast(case.m[,c(1,2,4)],genus~type,sum,fill = 0)
control.m2 = dcast(control.m[,c(1,2,4)],genus~type,sum,fill = 0)
gg = intersect(as.character(case.m2$genus),as.character(control.m2$genus))

ob1 = ob[as.character(ob$taxon_names) %in% as.character(gg),c(1,12)]
cr1 = cr[as.character(cr$taxon_names) %in% as.character(gg),c(1,12)]
crc1 = crc[as.character(crc$taxon_names) %in% as.character(gg),c(1,12)]
asd1 = asd[as.character(asd$taxon_names) %in% as.character(gg),c(1,12)]
aids1 = aids[as.character(aids$taxon_names) %in% as.character(gg),c(1,12)]

ob1$type = "OB"
cr1$type = "CR"
crc1$type = "CRC"
asd1$type = "ASD"
aids1$type = "AIDS"

mydata = rbind(ob1,asd1,aids1,cr1,crc1)
mydata[which(is.na(mydata$netzgo.scale)),"netzgo.scale"] = 0
mydata2 = dcast(mydata[,c(1,3,2)],taxon_names~type,sum,fill = 0)
mydata.single = mydata2


##all
mydata.multi2 = melt(mydata.multi,id.vars = "taxon_names")
mydata.multi2$type = "multiple"
mydata.single2 = melt(mydata.single,id.vars = "taxon_names")
mydata.single2$type = "single"
mydata.all = rbind(mydata.multi2,mydata.single2)

dis = unique(mydata.all$variable)
a1 = mydata.all[which(mydata.all$type == "multiple"),]
a2 = mydata.all[which(mydata.all$type == "single"),]
m.data = data.frame()
for (n in 1:length(dis))
{
  m.data[n,1] = dis[n]
  m.data[n,2] = "multiple"
  m.data[n,3] = mean(as.numeric(as.character(a1[which(a1$variable == as.character(dis[n])),"value"])))
  m.data[n,4] = sd(as.numeric(as.character(a1[which(a1$variable == as.character(dis[n])),"value"])))
}
colnames(m.data) = c("disease","type","mean","sd")  
s.data = data.frame()
for (n in 1:length(dis))
{
  s.data[n,1] = dis[n]
  s.data[n,2] = "single"
  s.data[n,3] = mean(as.numeric(as.character(a2[which(a2$variable == as.character(dis[n])),"value"])))
  s.data[n,4] = sd(as.numeric(as.character(a2[which(a2$variable == as.character(dis[n])),"value"])))
}
colnames(s.data) = c("disease","type","mean","sd")  

####plot fig5e
a3 = rbind(m.data,s.data)
a4 = dcast(a3[,1:3],type~disease,mean,fill = 0)
ggradar(a4,grid.max = 0.5,grid.min = 0,values.radar = c("0","0.5"))

```

