Construct cross-species networkΒΆ

3.1_CrossSpeciesNetwork.v2.utf8
library(tidyverse)
library(igraph)
library(openxlsx)
library(stats)

Load files

t = 0.2
load("2.1.c_WGCNA_ThresholdingPower.UT.RData")
load("2.2.c_WGCNA.Modules.UT.RData")
load("2.1.z_WGCNA_ThresholdingPower.UT.RData")
load("2.2.z_WGCNA.Modules.UT.RData")

Construction of cross-species network

rename = function(mat){
  x2 = rownames(mat)
  x2[rownames(mat) == "201602_Is1"] = "201602_Mi1"
  x2[rownames(mat) == "201602_Mi1"] = "201602_Is1"
  rownames(mat) = x2
  mat = mat[sort(rownames(mat)), ]
  mat
}
eig.c = modules.c[[as.character(t)]]$newMEs
eig.z = modules.z[[as.character(t)]]$newMEs
eig.c = rename(eig.c)
eig.z = rename(eig.z)

colnames(eig.c) = paste0(gsub("^ME", "", colnames(eig.c)), "_c")
colnames(eig.z) = paste0(gsub("^ME", "", colnames(eig.z)), "_z")

cormat.cc = cor(eig.c, eig.c, method="pearson")
cormat.zz = cor(eig.z, eig.z, method="pearson")
cormat.cz = cor(eig.c, eig.z, method="pearson")

convert_cormat = function(mat){
  network = matrix(nrow=0, ncol=5)
  colnames(network) = c("Node1", "Node2", "Correlation", "Edge.color", "Edge.width")
  for(i in 1:(nrow(mat)-1)){
    for(j in (i+1):ncol(mat)){
      val = mat[i, j]
      if(abs(val) > 0.7){
        network = rbind(network, 
                        c(rownames(mat)[i], colnames(mat)[j], val, 
                          ifelse(val > 0, "#F46D43", "#74ADD1"),
                          ifelse(abs(val) > 0.85, 3, 1)))
      }
    }
  }
  network = data.frame(network, stringsAsFactors=F)
  network$Correlation = as.numeric(network$Correlation)
  network$Edge.width = as.numeric(network$Edge.width)
  network
}

network = rbind(convert_cormat(cormat.cc), 
                convert_cormat(cormat.zz))

tmp = cormat.cz %>% as.data.frame() %>% 
  rownames_to_column("Node1") %>%
  gather(key="Node2", value="Correlation", -Node1) %>% 
  mutate(Edge.color = ifelse(Correlation > 0, "#F46D43", "#74ADD1"),
         Edge.width = ifelse(abs(Correlation) > 0.9, 7,
                             ifelse(abs(Correlation) > 0.7, 5, 
                                    ifelse(Correlation > 0.5, 3, 1)))) %>%
  filter(abs(Correlation) > 0.3)
network = rbind(network, tmp)
network = network %>% mutate(Abs_Correlation = abs(Correlation), Pos=ifelse(Correlation > 0, 1, 0))
network %>% arrange(desc(Correlation)) %>% head()
##           Node1          Node2 Correlation Edge.color Edge.width
## 1    skyblue3_c          red_z   0.9418489    #F46D43          7
## 2      brown4_c     thistle2_z   0.8771784    #F46D43          5
## 3 saddlebrown_c    royalblue_z   0.8673096    #F46D43          5
## 4        cyan_c        plum2_z   0.8547272    #F46D43          5
## 5         tan_c         cyan_z   0.8534268    #F46D43          5
## 6 darkorange2_c navajowhite2_z   0.8397811    #F46D43          5
##   Abs_Correlation Pos
## 1       0.9418489   1
## 2       0.8771784   1
## 3       0.8673096   1
## 4       0.8547272   1
## 5       0.8534268   1
## 6       0.8397811   1

Correlation with zooxanthellae density

load("../Enviroment/ZooxanthellaeDensity.RData")
matZooDen = matZooDen %>% mutate(ID = paste0(str_extract(ID, "^[0-9]+"), "_", str_extract(ID, "...$")))

tmp = eig.c %>% rownames_to_column("ID") %>% left_join(matZooDen, by="ID") %>% column_to_rownames("ID") %>% na.omit()
cor.zoo.c = sapply(colnames(eig.c), 
                   function(x){ 
                     a = cor.test(tmp[[x]], tmp$Value, method="pearson")
                     c(a$estimate, pvalue=a$p.value)
                     })
cor.zoo.c = t(cor.zoo.c) %>% as.data.frame() %>% rownames_to_column("modules") %>% 
  mutate(FDR = p.adjust(pvalue, method="BH")) %>% filter(FDR < 0.1) %>% arrange(desc(cor))
cor.zoo.c
##          modules        cor      pvalue        FDR
## 1   orangered4_c  0.4222287 0.002026715 0.04219109
## 2      sienna3_c  0.4168325 0.002343949 0.04219109
## 3   lightcyan1_c  0.4009956 0.003543787 0.04252545
## 4     skyblue3_c  0.3674112 0.007995450 0.05756724
## 5         blue_c  0.3468131 0.012657735 0.06863775
## 6      bisque4_c  0.3318809 0.017349617 0.06863775
## 7    royalblue_c  0.3286855 0.018525931 0.06863775
## 8  saddlebrown_c  0.3272757 0.019066042 0.06863775
## 9       brown4_c -0.3373235 0.015492079 0.06863775
## 10     skyblue_c -0.3677189 0.007939003 0.05756724
tmp = eig.z %>% rownames_to_column("ID") %>% left_join(matZooDen, by="ID") %>% column_to_rownames("ID") %>% na.omit()
cor.zoo.z = sapply(colnames(eig.z), 
                   function(x){ 
                     a = cor.test(tmp[[x]], tmp$Value, method="pearson")
                     c(a$estimate, pvalue=a$p.value)
                     })
cor.zoo.z = t(cor.zoo.z) %>% as.data.frame() %>% rownames_to_column("modules") %>% 
  mutate(FDR = p.adjust(pvalue, method="BH")) %>% filter(FDR < 0.1) %>% arrange(desc(cor))
cor.zoo.z
##             modules        cor       pvalue          FDR
## 1   darkturquoise_z  0.5874042 5.860913e-06 0.0001641056
## 2          brown4_z  0.5135533 1.161720e-04 0.0010842722
## 3           green_z  0.3268130 1.924618e-02 0.0673616164
## 4             red_z  0.3126611 2.549468e-02 0.0793167811
## 5        thistle2_z -0.3512329 1.149742e-02 0.0459896865
## 6         salmon4_z -0.3594602 9.579661e-03 0.0447050858
## 7        skyblue3_z -0.4359049 1.386890e-03 0.0077665857
## 8           plum2_z -0.4777099 3.926541e-04 0.0027485786
## 9 lightsteelblue1_z -0.5390238 4.485743e-05 0.0006280041

Export network as Cytoscape format files

# Save for Cytoscape
modules = c(cor.zoo.c$modules, cor.zoo.z$modules)
network = network %>% filter(Node1 %in% modules) %>% filter(Node2 %in% modules)
write.table(network, file="3.1_CrossSpeciesNetwork.Edge.v2.txt", sep="\t", row.names=F, quote=F)

# Save edge information
node.c = data.frame(color = modules.c[[as.character(t)]]$colors, stringsAsFactors=F) %>% 
  group_by(color) %>% count(name="size") %>%
  mutate(modules = paste0(color, "_c"), species="coral", hexcol=gplots::col2hex(color)) %>%
  left_join(cor.zoo.c, by="modules") %>% na.omit() %>%
  mutate(group = ifelse(cor > 0, "coral_1", "coral_2")) %>% 
  select(color, size, modules, species, hexcol, group)
node.z = data.frame(color = modules.z[[as.character(t)]]$colors, stringsAsFactors=F) %>% 
  group_by(color) %>% count(name="size") %>%
  mutate(modules = paste0(color, "_z"), species="zooxanthellae", hexcol=gplots::col2hex(color)) %>%
  left_join(cor.zoo.z, by="modules") %>% na.omit() %>%
  mutate(group = ifelse(cor > 0, "zoo_1", "zoo_2")) %>% 
  select(color, size, modules, species, hexcol, group)
node = rbind(node.c, node.z)
write.table(node, file="3.1_CrossSpeciesNetwork.Node.v2.txt", sep="\t", row.names=F, quote=F)