7. co-occurrance analysis.Rmd

---
title: "7. co-occurrance networks"
author: "Baptiste Oosterlinck"
date: "3/30/2022"
output: html_document
---

#function to extract and format data from the phyloseq object for the cooccur analysis
```{r}
NetworkPrepare <- function(ps){
  #getting the relative abundance table for the tumor samples
  abundance <- otu_table(ps) %>%
    as.data.frame() %>%
    mutate(.,ASV = row.names(.))
  #replace the ASV-names by the genus names
  taxon <- tax_table(ps) %>%
    as.data.frame() %>%
    mutate(.,ASV = row.names(.))%>%
    dplyr::select(.,c("ASV","Genus"))
  #merging the abundance and taxonomy into one
  abundance <- left_join(abundance,taxon,by = "ASV",keep = FALSE) %>%
    mutate(., ASV = NULL)%>%
    column_to_rownames(.,var = "Genus")
  
  abundance[abundance > 0] <- 1
  return(abundance)
} 
```

##function to extract and visualise the co-occurance results
```{r}
NetworkSig <- function(Network){
  #extracting significant results
  Network_Sig <- Network[["results"]] %>%
    subset(., p_lt < 0.05 | p_gt < 0.05)
  #getting the names and IDs
  Node_Names <- unique(c(Network_Sig$sp1_name,Network_Sig$sp2_name))
  Node_IDS <- unique(c(Network_Sig$sp1,Network_Sig$sp2))
  #creating the nodes dataframe
  nodes <- data.frame(id= Node_IDS, 
                    label = Node_Names,
                    color = "#606482",
                    shadow = TRUE)
  #creating the edges dataframe
  edges <- data.frame(from = Network_Sig$sp1, to = Network_Sig$sp2,
                    color = ifelse(Network_Sig$p_lt <= 0.05, "#eba417","#069910"),
                    dashes = ifelse(Network_Sig$p_lt <= 0.05,yes = TRUE, no = FALSE))
  return(list(Significant = Network_Sig,nodes = nodes, edges = edges))
}
```

#1. Building the dataframes to do the co-occurance analysis
##1.1 doing the co-occurance per tissue type control - tumor
```{r}
#getting the presence absence data from the tumor samples
adundance_T <- subset_samples(agglom_relAb2, Tissue.Type == "tumor") %>%
  tax_glom(.,"Genus") %>%
  t(.)%>%
  NetworkPrepare(.)

#running the co-occurrence for the tumor samples
Network_Tumor <- cooccur(mat = adundance_T,type = "spp_site",spp_names = TRUE)

#getting the presence absence data from the control samples
abundance_C <-subset_samples(agglom_relAb2, Tissue.Type == "noninflammed") %>%
  tax_glom(.,"Genus") %>%
  t(.) %>%
  NetworkPrepare(.)

#running the co-occurrence for the control samples
Network_Control <- cooccur(mat = abundance_C,type = "spp_site",spp_names = TRUE)
```

##1.2 creating the visualistation of the co-occurance in the control samples
```{r}
#extracting the results from the co-occurrance object
  #filtering the significant results from the dataframe:
NetCont_resSig <- Network_Control[["results"]] %>%
  mutate(.,cohort = "control") %>%
  subset(., p_lt < 0.05 | p_gt < 0.05)

Node_Names <- unique(c(NetCont_resSig$sp1_name,NetCont_resSig$sp2_name))
Node_IDS <- unique(c(NetCont_resSig$sp1,NetCont_resSig$sp2))
#creating the nodes of the network
nodes <- data.frame(id= Node_IDS, 
                    label = Node_Names,
                    color = "#606482",
                    shadow = TRUE)
head(nodes)

edges <- data.frame(from = NetCont_resSig$sp1, to = NetCont_resSig$sp2,
                    color = ifelse(NetCont_resSig$p_lt <= 0.05, "#eba417","#069910"),
                    dashes = ifelse(NetCont_resSig$p_lt <= 0.05,TRUE,FALSE))
summary(edges)

head(edges)
#plotting the network
visNetwork(nodes = nodes, edges = edges) %>%
  visIgraphLayout(layout = "layout_nicely")
```

##1.3 creating the visualistation of the co-occurance in the Tumor samples
```{r}
NetTum_resSig <- Network_Tumor[["results"]] %>%
  mutate(.,cohort = "tumor")%>%
  subset(., p_lt < 0.05 | p_gt < 0.05)

Node_Names <- unique(c(NetTum_resSig$sp1_name,NetTum_resSig$sp2_name))
Node_IDS <- unique(c(NetTum_resSig$sp1,NetTum_resSig$sp2))
#creating the nodes of the network
nodes <- data.frame(id= Node_IDS, 
                    label = Node_Names,
                    color = "#606482",
                    shadow = TRUE)
head(nodes)

edges <- data.frame(from = NetTum_resSig$sp1, to = NetTum_resSig$sp2,
                    color = ifelse(NetTum_resSig$p_lt <= 0.05, "#eba417","#069910"),
                    dashes = ifelse(NetTum_resSig$p_lt <= 0.05,TRUE,FALSE))
summary(edges)

head(edges)
#plotting the network
visNetwork(nodes = nodes, edges = edges) %>%
  visIgraphLayout(layout = "layout_with_kk")
```

#2. Getting co-occurance per mucin phenotype
```{r}
###########Gastric phenotye#############################
#getting presence absance data from the tumor agglom_BAT
adundance_T_Gastric <- subset_samples(agglom_BAT, Mucin.Phenotype == "Gastric") %>%
  t(.) %>%
  NetworkPrepare(.)
#running the cooccurance  
Network_T_Gastric <- cooccur(mat = adundance_T_Gastric,type = "spp_site",spp_names = TRUE)

###########Intestinal phenotye#############################
#getting presence absance data from the tumor agglom_BAT
adundance_T_Intestinal <- subset_samples(agglom_BAT, Mucin.Phenotype == "Intestinal") %>%
  t(.) %>%
  NetworkPrepare(.)
#running the cooccurance  
Network_T_Intestinal <- cooccur(mat = adundance_T_Intestinal,type = "spp_site",spp_names = TRUE)

###########Mixed phenotye#############################
#getting presence absance data from the tumor agglom_BAT
adundance_T_Mixed <- subset_samples(agglom_BAT, Mucin.Phenotype == "Mixed") %>%
  t(.) %>%
  NetworkPrepare(.)
#running the cooccurance  
Network_T_Mixed <- cooccur(mat = adundance_T_Mixed,type = "spp_site",spp_names = TRUE)

###########Null phenotye#############################
#getting presence absance data from the tumor agglom_BAT
adundance_T_Null <- subset_samples(agglom_BAT, Mucin.Phenotype == "Null") %>%
  t(.) %>%
  NetworkPrepare(.)
#running the cooccurance  
Network_T_Null <- cooccur(mat = adundance_T_Null,type = "spp_site",spp_names = TRUE)
```

##2.1 visualising the networks
```{r}
#getting all differential abundant genera from the aldex objects:
Differential_genera <- data.frame(Genus = AldexMUC5AC_Tax@.Data[,"Genus"],Mucin = "MUC5AC") %>%
  rbind(., data.frame(Genus = AldexMUC6_Tax@.Data[,"Genus"],Mucin = "MUC6")) %>%
  rbind(.,data.frame(Genus = AldexMUC2_Tax@.Data[,"Genus"],Mucin = "MUC2")) %>%
  rbind(.,data.frame(Genus = AldexMUC4_Tax@.Data[,"Genus"],Mucin = "MUC4")) %>%
  rbind(.,data.frame(Genus = AldexMUC13_Tax@.Data[,"Genus"],Mucin = "MUC13"))
Differential_genera$EnrichedGroup <- c("Low","High","High","High","High","High","Mid","Mid","Mid","High","Low","High","High","High","Low","High","Low","High","Low","Low","High","Low")


#function to change the node shape according to the mucin it associates with
NodeShape <- function(x, Diff_gen){
  shape <- c()
  for(label in x$nodes$label){
  if(label %in% Diff_gen$Genus){
    if(Diff_gen[Diff_gen$Genus == label, "Mucin"] == "MUC13"){
      shape <- c(shape, "MUC13")
    }else if(Diff_gen[Diff_gen$Genus == label, "Mucin"] == "MUC4"){
      shape <- c(shape, "MUC4")
    }else if(Diff_gen[Diff_gen$Genus == label, "Mucin"] == "MUC2"){
      shape <- c(shape, "MUC2")
    }else if(Diff_gen[Diff_gen$Genus == label, "Mucin"] == "MUC6"){
      shape <- c(shape,"MUC6")
    }else if(Diff_gen[Diff_gen$Genus == label, "Mucin"] == "MUC5AC"){
      shape <- c(shape,"MUC5AC")
    }else{
      shape <- c(shape,"NA")
    }
  }else{
    shape <- c(shape,"NA")
  }
  }
  return(shape)
}

#function to change the node color according to the expression level it associates with
NodeColor <- function(x, Diff_gen){
  color <- c()
  for(label in x$nodes$label){
    if(label %in% Diff_gen$Genus){
      if(Diff_gen[Diff_gen$Genus == label,"EnrichedGroup"] == "High"){
        color <- c(color, "High")
      }else if(Diff_gen[Diff_gen$Genus == label,"EnrichedGroup"] == "Mid"){
        color <- c(color,"Mid")
      }else if(Diff_gen[Diff_gen$Genus == label,"EnrichedGroup"] == "Low"){
        color <- c(color,"Low")
      }
    }else
      color <- c(color,"Not differentially abundant")
  }
  return(color)
}

##preparing the network of the gastric mucin phenotype for visualization
Gastric <- NetworkSig(Network_T_Gastric)
Differential_Gastric <- subset(Differential_genera, Mucin %in% c("MUC5AC","MUC6"))
Gastric$nodes$shape <- NodeShape(Gastric, Differential_Gastric)
Gastric$nodes$color <- NodeColor(Gastric, Differential_Gastric)

#using network package
Gastric_network <- network(Gastric$edges, vertex.attr = Gastric$nodes, matrix.type = "edgelist",ignore.eval = FALSE,directed = FALSE) %>%
  as_tbl_graph
#plotting the network
Gastric_NetworkGraph <- ggraph(Gastric_network,layout = "linear",circular = TRUE)+
  geom_edge_link(aes(linetype = dashes, colour = dashes))+
  geom_node_point(aes(shape = shape, colour = color), size = 3)+
  geom_node_text(aes(label = label, colour = color), size = 1,repel = TRUE, max.overlaps = 50)+
  theme_graph(plot_margin = margin(0.02,0.02,0.02,0.02,"cm"))+
  scale_color_manual(limits = c("High","Mid","Low"),values = c("#619CFF","#00BA38","#F8766D"))+
  scale_shape_manual(limits = c("MUC13","MUC2","MUC4","MUC5AC","MUC6","NA"), values = c(2,5,0,10,12,20))+
  ggtitle("A) Network of gastric mucin phenotype")

##preparing the network of the intestinal mucin phenotype for visualization
Intestinal <- NetworkSig(Network_T_Intestinal)
Differential_Intestinal <- subset(Differential_genera, Mucin %in% c("MUC2","MUC4","MUC13"))
Intestinal$nodes$shape <- NodeShape(Intestinal, Differential_Intestinal)
Intestinal$nodes$color <- NodeColor(Intestinal, Differential_Intestinal)

#making network:
Intestinal_network <- network(Intestinal$edges, vertex.attr = Intestinal$nodes, matrix.type = "edgelist",ignore.eval = FALSE, directed = FALSE) %>%
  as_tbl_graph

#plotting the network
Intestinal_NetworkGraph <- ggraph(Intestinal_network,layout = "linear", circular=TRUE)+
  geom_edge_link(aes(linetype = dashes,colour = dashes))+
  geom_node_point(aes(shape = shape, colour = color), size = 3)+
  geom_node_text(aes(label = label, colour = color), size = 1,repel = TRUE, max.overlaps = 50)+
  theme_graph(plot_margin = margin(0.02,0.02,0.02,0.02,"cm"))+
  scale_color_manual(limits = c("High","Mid","Low"),values = c("#619CFF","#00BA38","#F8766D"))+
  scale_shape_manual(limits = c("MUC13","MUC2","MUC4","MUC5AC","MUC6","NA"), values = c(2,5,0,10,12,20))+
  ggtitle("B) Network of intestinal mucin phenotype")

##preparing the network of the Mixed mucin phenotype for visualization
Mixed <- NetworkSig(Network_T_Mixed)
Mixed$nodes$shape <- NodeShape(Mixed, Differential_genera)
Mixed$nodes$color <- NodeColor(Mixed, Differential_genera)
Mixed$nodes$size <- NodeSize(Mixed, Differential_genera)

Mixed_network <- network(Mixed$edges, vertex.attr = Mixed$nodes, matrix.type = "edgelist",ignore.eval = FALSE, directed = FALSE) %>%
  as_tbl_graph

#plotting the network
Mixed_NetworkGraph <- ggraph(Mixed_network,layout = "linear", circular=TRUE)+
  geom_edge_link(aes(linetype = dashes,colour = dashes))+
  geom_node_point(aes(shape = shape, colour = color),size = 3)+
  geom_node_text(aes(label = label, colour = color), size = 1,repel = TRUE, max.overlaps = 50)+
  theme_graph(plot_margin = margin(0.02,0.02,0.02,0.02,"cm"))+
  scale_color_manual(limits = c("High","Mid","Low"),values = c("#619CFF","#00BA38","#F8766D"))+
  scale_shape_manual(limits = c("MUC13","MUC2","MUC4","MUC5AC","MUC6","NA"), values = c(2,5,0,10,12,20))+
  ggtitle("C) Network of mixed mucin phenotype")

##preparing the network of the null mucin phenotype for visualization
Null <- NetworkSig(Network_T_Null)
Null$nodes$shape <- NodeShape(Null, Differential_genera)
Null$nodes$color <- NodeColor(Null, Differential_genera)
Null$nodes$size <- NodeSize(Null, Differential_genera)

Null_network <- network(Null$edges, vertex.attr = Null$nodes, matrix.type = "edgelist",ignore.eval = FALSE, directed = FALSE) %>%
  as_tbl_graph

Null_NetworkGraph <- ggraph(Null_network,layout = "linear", circular=TRUE)+
  geom_edge_link(aes(linetype = dashes,colour = dashes))+
  geom_node_point(aes(shape = shape, colour = color),size = 3)+
  geom_node_text(aes(label = label, colour = color),size = 1,repel = TRUE,max.overlaps = 50)+
  theme_graph(plot_margin = margin(0.02,0.02,0.02,0.02,"cm"))+
  scale_color_manual(limits = c("High","Mid","Low"),values = c("#619CFF","#00BA38","#F8766D"))+
  scale_shape_manual(limits = c("MUC13","MUC2","MUC4","MUC5AC","MUC6","NA"), values = c(2,5,0,10,12,20))+
  ggtitle("D) Network of null mucin phenotype")

ggarrange(Gastric_NetworkGraph,Intestinal_NetworkGraph,Mixed_NetworkGraph,Null_NetworkGraph, ncol = 2,nrow = 2, common.legend = TRUE,legend = "right")
```

##2.2 quantification of the network interactions
```{r}
#total number of significant interactions
proportionSignif <- function(network){
  analyzed_temp <- network[["pairs"]]
  NetSig_temp <- NetworkSig(network)$Significant
  analyzedSig_temp <- nrow(NetSig_temp)
  pos_temp <- nrow(subset(NetSig_temp, p_gt <= 0.05))
  neg_temp <- nrow(subset(NetSig_temp, p_lt <= 0.05))

  data.frame(general = analyzedSig_temp, possitive = pos_temp, negative = neg_temp)
}

Network_List <- list(Network_T_Gastric, Network_T_Intestinal,Network_T_Mixed, Network_T_Null)
names(Network_List) <- c("Gastric","Intestinal","Mixed","Null")

prop_Signif <- lapply(Network_List, proportionSignif)
names(prop_Signif) <- c("Gastric","Intestinal","Mixed","Null")
prop_Signif <- do.call(rbind.data.frame,prop_Signif) %>%
  transpose()
rownames(prop_Signif) <- c("general", "positive","negative")
colnames(prop_Signif) <- c("Gastric","Intestinal","Mixed","Null")

PropSignif_Chi2 <- apply(prop_Signif,MARGIN = 1, chisq.test)

#visualising the Pearson residuals of the chisquare test
residuals <- PropSignif_Chi2$residuals%>%
  as.data.frame %>%
  rownames_to_column(var = "connection_type") %>%
  gather(., key = Phenotype, value = residual,2:5)
ggplot(data = residuals, aes(x= Phenotype, y = connection_type, fill = residual))+
  geom_tile()

#difference in number of interactions of a genus per phenotype
NetworkProp <- function(network){
  networkSig_temp <- NetworkSig(network)
  Bacteria_Temp <- unique(c(networkSig_temp$Significant$sp1_name,networkSig_temp$Significant$sp2_name))
  total_connections <- nrow(networkSig_temp$Significant)
  Overview_connections <- data.frame(genus = c(), proportion = c())
  for(bact in Bacteria_Temp){
    bact_connections <- nrow(subset(networkSig_temp$Significant, sp1_name == bact))
    bact_poss <- nrow(subset(networkSig_temp$Significant, sp1_name == bact & p_gt <= 0.05))
    bact_neg <- nrow(subset(networkSig_temp$Significant, sp1_name == bact & p_lt <= 0.05))
    Overview_connections <- rbind(Overview_connections, c(bact, bact_connections, bact_poss, bact_neg))
  }
  colnames(Overview_connections) <- c("Genus","General_interactions","Pos_interactions","Neg_interactions")
  Overview_connections
  }

Proportions_Interactions <- lapply(Network_List, NetworkProp)
Proportions_Interactions <- reduce(Proportions_Interactions, full_join, by = "Genus")
colnames(Proportions_Interactions) <- c('Genus',"Gastric_G","Gastric_P","Gastric_N","Intestinal_G","Intestinal_P","Intestinal_N","Mixed_G","Mixed_P","Mixed_N","Null_G","Null_P","Null_N")
Proportions_Interactions[,2:13] <- apply(Proportions_Interactions[,2:13],MARGIN = 2, as.numeric)

#replace NA by 0
Proportions_Interactions[is.na(Proportions_Interactions)] <- 0
#remove all rows that sum up to 0
Proportions_Interactions <- subset(Proportions_Interactions,rowSums(Proportions_Interactions[,2:13]) != 0)
rownames(Proportions_Interactions) <- Proportions_Interactions$Genus
Proportions_Interactions$Genus <- NULL

GenInteractions_Chi2 <- chisq.test(Proportions_Interactions[,c("Gastric_G","Intestinal_G","Mixed_G","Null_G")])

Proportions_Interactions <- GenInteractions_Chi2$residuals %>%
  as.data.frame %>%
  rownames_to_column(var = "Genus") %>%
  full_join(rownames_to_column(Proportions_Interactions,var = "Genus"),.,by = "Genus")

PosInteractions_Chi2 <- chisq.test(Proportions_Interactions[,c("Gastric_P","Intestinal_P","Mixed_P","Null_P")]%>%subset(.,rowSums(.) != 0))

Proportions_Interactions <- PosInteractions_Chi2$residuals %>%
  as.data.frame %>%
  rownames_to_column(var = "Genus") %>%
  full_join(Proportions_Interactions,.,by = "Genus")

NegInteractions_Chi2 <- chisq.test(Proportions_Interactions[,c("Gastric_N","Intestinal_N","Mixed_N","Null_N")]%>%subset(.,rowSums(.) != 0))

Proportions_Interactions <- NegInteractions_Chi2$residuals %>%
  as.data.frame %>%
  rownames_to_column(var = "Genus") %>%
  full_join(Proportions_Interactions,.,by = "Genus")

writexl::write_xlsx(Proportions_Interactions, path = "Phenotype_networks_chisq.xlsx")
```

#3. Getting cooccurance for MUC13 high - low

```{r}
###########Low MUC13#############################
#getting presence/absence data from the tumor agglom_BAT
adundance_T_MUC13Low <- subset_samples(agglom_BAT, Strat.Score_MUC13_T == -1) %>%
  t(.) %>%
  NetworkPrepare(.)

#running the cooccurance  
Network_MUC13Low <- cooccur(mat =adundance_T_MUC13Low,type = "spp_site",spp_names = TRUE)

MUC13Low <- NetworkSig(Network_MUC13Low)
Differential_MUC13 <- subset(Differential_genera,Mucin == "MUC13")
MUC13Low$nodes$shape <- NodeShape(MUC13Low, Differential_MUC13)
MUC13Low$nodes$color <- NodeColor(MUC13Low, Differential_MUC13)
MUC13Low$nodes$size <- NodeSize(MUC13Low, Differential_MUC13)

MUC13Low_Network <- network(MUC13Low$edges, vertex.attr = MUC13Low$nodes, matrix.type = "edgelist",ignore.eval = FALSE, directed = FALSE) %>%
  as_tbl_graph

MUC13Low_NetworkGraph <- ggraph(MUC13Low_Network,layout = "linear", circular=TRUE)+
  geom_edge_link(aes(linetype = dashes,colour = dashes))+
  geom_node_point(aes(shape = shape, colour = color),size = 3)+
  geom_node_text(aes(label = label, colour = color),size = 1,repel = TRUE,max.overlaps = 50)+
  theme_graph(plot_margin = margin(0.02,0.02,0.02,0.02,"cm"))+
  scale_color_manual(limits = c("High","Mid","Low"),values = c("#619CFF","#00BA38","#F8766D"))+
  scale_shape_manual(limits = c("MUC13","MUC2","MUC4","MUC5AC","MUC6","NA"), values = c(2,5,0,10,12,20))+
  ggtitle("A) Network of Low MUC13 tumor tissues")

###########Hihg MUC13#############################
#getting presence absance data from the tumor agglom_BAT
adundance_T_MUC13High <- subset_samples(agglom_BAT, Strat.Score_MUC13_T == 1) %>%
  t(.) %>%
  NetworkPrepare(.)

#running the cooccurance  
Network_MUC13High <- cooccur(mat =adundance_T_MUC13High,type = "spp_site",spp_names = TRUE)

MUC13High <- NetworkSig(Network_MUC13High)
MUC13High$nodes$shape <- NodeShape(MUC13High, Differential_MUC13)
MUC13High$nodes$color <- NodeColor(MUC13High, Differential_MUC13)
MUC13High$nodes$size <- NodeSize(MUC13High, Differential_MUC13)

MUC13High_Network <- network(MUC13High$edges, vertex.attr = MUC13High$nodes, matrix.type = "edgelist",ignore.eval = FALSE, directed = FALSE) %>%
  as_tbl_graph

MUC13High_NetworkGraph <- ggraph(MUC13High_Network,layout = "linear", circular=TRUE)+
  geom_edge_link(aes(linetype = dashes,colour = dashes))+
  geom_node_point(aes(shape = shape, colour = color),size = 3)+
  geom_node_text(aes(label = label, colour = color),size = 1,repel = TRUE,max.overlaps = 50)+
  theme_graph(plot_margin = margin(0.02,0.02,0.02,0.02,"cm"))+
  scale_color_manual(limits = c("High","Mid","Low"),values = c("#619CFF","#00BA38","#F8766D"))+
  scale_shape_manual(limits = c("MUC13","MUC2","MUC4","MUC5AC","MUC6","NA"), values = c(2,5,0,10,12,20))+
  ggtitle("A) Network of High MUC13 tumor tissues")

ggarrange(MUC13Low_NetworkGraph,MUC13High_NetworkGraph,ncol = 2, common.legend = TRUE,legend = "right")
```

#3.1 quantification of the MUC13 networks
```{r}
Network_List <- list(Network_MUC13Low, Network_MUC13High)

prop_Signif <- lapply(Network_List, proportionSignif)
names(prop_Signif) <- c("MUC13_Low","MUC13_High")
prop_Signif <- do.call(rbind.data.frame,prop_Signif) %>%
  transpose()
rownames(prop_Signif) <- c("general", "positive","negative")
colnames(prop_Signif) <- c("MUC13_Low","MUC13_High")

PropSignif_Chi2 <- apply(prop_Signif,MARGIN = 1, chisq.test)

#visualising the Pearson residuals of the chisquare test
residuals <- PropSignif_Chi2$residuals%>%
  as.data.frame %>%
  rownames_to_column(var = "connection_type") %>%
  gather(., key = Phenotype, value = residual,2:3)
ggplot(data = residuals, aes(x= Phenotype, y = connection_type, fill = residual))+
  geom_tile()

#difference in number of interactions of a genus per phenotype
Proportions_Interactions <- lapply(Network_List, NetworkProp)
names(Proportions_Interactions) <- c("MUC13_Low","MUC13_High")
Proportions_Interactions <- reduce(Proportions_Interactions, full_join, by = "Genus")
colnames(Proportions_Interactions) <- c("Genus","MUC13_Low_G","MUC13_Low_P","MUC13_Low_N","MUC13_High_G","MUC13_High_P","MUC13_High_N")
Proportions_Interactions[,2:7] <- apply(Proportions_Interactions[,2:7],MARGIN = 2, as.numeric)

#replace NA by 0
Proportions_Interactions[is.na(Proportions_Interactions)] <- 0
rownames(Proportions_Interactions) <- Proportions_Interactions$Genus

GenInteractions_Chi2 <- chisq.test(Proportions_Interactions[,c("MUC13_Low_G","MUC13_High_G")]%>%
                                     subset(.,rowSums(.) != 0))

Proportions_Interactions <- GenInteractions_Chi2$residuals %>%
  as.data.frame %>%
  rownames_to_column(var = "Genus") %>%
  full_join(Proportions_Interactions,.,by = "Genus")

rownames(Proportions_Interactions) <- Proportions_Interactions$Genus
PosInteractions_Chi2 <- chisq.test(Proportions_Interactions[,c("MUC13_Low_P","MUC13_High_P")]
                                   %>%subset(.,rowSums(.) != 0))

Proportions_Interactions <- PosInteractions_Chi2$residuals %>%
  as.data.frame %>%
  rownames_to_column(var = "Genus") %>%
  full_join(Proportions_Interactions,.,by = "Genus")

rownames(Proportions_Interactions) <- Proportions_Interactions$Genus
NegInteractions_Chi2 <- chisq.test(Proportions_Interactions[,c("MUC13_Low_N","MUC13_High_N")]%>%
                                     subset(.,rowSums(.) != 0))

Proportions_Interactions <- NegInteractions_Chi2$residuals %>%
  as.data.frame %>%
  rownames_to_column(var = "Genus") %>%
  full_join(Proportions_Interactions,.,by = "Genus")

writexl::write_xlsx(Proportions_Interactions, path = "tabellen/MUC13_networks_chisq.xlsx")
```