---
Author: "Amin Haghani"
Paper: "DNA Methylation Networks Underlying Mammalian Traits"
Authors: "Amin Haghani1, 2 †*; Caesar Z. Li3, 4 †; Todd R. Robeck5; Joshua Zhang1; Ake T. Lu1, 2; Julia Ablaeva6; Victoria A. Acosta-Rodríguez7; Danielle M. Adams8; Abdulaziz N. Alagaili9, 10; Javier Almunia11; Ajoy Aloysius12; Nabil M.S. Amor13; Reza Ardehali14; Adriana Arneson15, 16; C. Scott Baker17; Gareth Banks18; Katherine Belov19; Nigel C. Bennett20; Peter Black21; Daniel T. Blumstein22, 23; Eleanor K. Bors17; Charles E. Breeze24; Robert T. Brooke25; Janine L. Brown26; Gerald Carter27; Alex Caulton28, 29; Julie M. Cavin30; Lisa Chakrabarti31; Ioulia Chatzistamou32; Andreas S. Chavez27, 33; Hao Chen34; Kaiyang Cheng35; Priscila Chiavellini36; Oi-Wa Choi37, 38; Shannon Clarke28; Joseph A. Cook39; Lisa N. Cooper40; Marie-Laurence Cossette41; Joanna Day42; Joseph DeYoung37, 38; Stacy Dirocco43; Christopher Dold44; Jonathan L. Dunnum39; Erin E. Ehmke45; Candice K. Emmons46; Stephan Emmrich6; Ebru Erbay47, 48, 49; Claire Erlacher-Reid43; Chris G. Faulkes50, 51; Zhe Fei3, 52; Steven H. Ferguson53, 54; Carrie J. Finno55; Jennifer E. Flower56; Jean-Michel Gaillard57; Eva Garde58; Livia Gerber59, 60; Vadim N. Gladyshev61; Rodolfo G. Goya36; Matthew J Grant62; Carla B. Green7; M. Bradley Hanson46; Daniel W. Hart20; Martin Haulena63; Kelsey Herrick64; Andrew N. Hogan65; Carolyn J. Hogg19; Timothy A. Hore66; Taosheng Huang67; Juan Carlos Izpisua Belmonte2; Anna J. Jasinska37, 68, 69; Gareth Jones70; Eve Jourdain71; Olga Kashpur72; Harold Katcher73; Etsuko Katsumata74; Vimala Kaza75; Hippokratis Kiaris76; Michael S. Kobor77; Pawel Kordowitzki78; William R. Koski79; Michael Krützen60; Soo Bin Kwon16, 15; Brenda Larison22, 80; Sang-Goo Lee61; Marianne Lehmann36; Jean-François Lemaître57; Andrew J. Levine81; Xinmin Li82; Cun Li83, 84; Andrea R. Lim1; David T.S. Lin85; Dana M. Lindemann43; Schuyler W.  Liphardt86; Thomas J. Little87; Nicholas Macoretta6; Dewey Maddox88; Craig O. Matkin89; Julie A. Mattison90; Matthew McClure91; June Mergl92; Jennifer J. Meudt93; Gisele A. Montano5; Khyobeni Mozhui94; Jason Munshi-South95; William J. Murphy96, 97; Asieh Naderi76; Martina Nagy98; Pritika Narayan62; Peter W. Nathanielsz83, 84; Ngoc B. Nguyen14; Christof Niehrs99, 100; Batsaikhan Nyamsuren101; Justine K. O'Brien42; Perrie O'Tierney Ginn72; Duncan T Odom102, 103; Alexander G. Ophir104; Steve Osborn105; Elaine A. Ostrander65; Kim M. Parsons46; Kimberly C. Paul81; Amy B. Pedersen87; Matteo Pellegrini106; Katharina J. Peters60, 107; Jessica L. Petersen108; Darren W. Pietersen109; Gabriela M. Pinho22; Jocelyn Plassais65; Jesse R. Poganik61; Natalia A. Prado110, 26; Pradeep Reddy111, 2; Benjamin Rey57; Beate R. Ritz112, 113, 81; Jooke Robbins114; Magdalena Rodriguez115; Jennifer Russell105; Elena Rydkina6; Lindsay L. Sailer104; Adam B. Salmon116; Akshay Sanghavi73; Kyle M. Schachtschneider117, 118, 119; Dennis Schmitt120; Todd Schmitt64; Lars Schomacher99; Lawrence B. Schook117, 121; Karen E. Sears22; Ashley W. Seifert12; Aaron B.A. Shafer122; Anastasia V. Shindyapina61; Melanie Simmons45; Kavita Singh123; Ishani Sinha22; Jesse Slone67; Russel G. Snell62; Elham Soltanmohammadi76; Matthew L. Spangler108; Maria Spriggs21; Lydia Staggs43; Nancy Stedman21; Karen J. Steinman124; Donald T Stewart125; Victoria J. Sugrue66; Balazs Szladovits126; Joseph S. Takahashi7, 127; Masaki Takasugi6; Emma C. Teeling128; Michael J. Thompson106; Bill Van Bonn129; Sonja C. Vernes130, 131; Diego Villar132; Harry V. Vinters133; Ha Vu15, 16; Mary C. Wallingford72; Nan Wang37, 38; Gerald S. Wilkinson8; Robert W. Williams134; Qi Yan3, 2; Mingjia Yao3; Brent G. Young54; Bohan Zhang61; Zhihui Zhang6; Yang Zhao6; Peng Zhao14, 135; Wanding Zhou136, 137; Joseph A. Zoller3; Jason Ernst15, 16; Andrei Seluanov138; Vera Gorbunova138; X. William Yang37, 38; Ken Raj139; Steve Horvath1, 2 *"
Date: "07-24-2023"
output: html_document
---

```{r libraries, message=FALSE}
library(easypackages)
libraries("readxl", "psych", "tidyr", "dplyr","gplots", "RColorBrewer","limma", "ggplot2", "metap", "anRichment", "clusterProfiler", "RColorBrewer", "ggstatsplot", "gridExtra", "ggupset", "psych", "corrplot", "limma", "edgeR", "gridExtra", "ggupset", "ggrepel", "data.table", "stringr")
```

```{r Network 1, reference}
net1MEs <- readRDS("WGCNA results/No marsupials/Module MEs merged net 1.RDS")
net1KMEs <- readRDS("WGCNA results/No marsupials/Merged modules KME net 1.RDS")
net1geneTree <- readRDS("WGCNA results/No marsupials/geneTree net 1.RDS")
net1mergedColors <- readRDS("WGCNA results/No marsupials/mergedColors net 1.RDS")
```

```{r Network 2}
net2MEs <- readRDS("WGCNA results/All, mappability filter/Module MEs merged net 2.RDS")
net2KMEs <- readRDS("WGCNA results/All, mappability filter/Merged modules KME net 2.RDS")
net2geneTree <- readRDS("WGCNA results/All, mappability filter/geneTree net 2.RDS")
net2mergedColors <- readRDS("WGCNA results/All, mappability filter/mergedColors net 2.RDS")


mergedData <- net1KMEs %>% tibble::rownames_to_column(var = "CGid") %>% dplyr::select(CGid) %>% mutate(net1 = net1mergedColors)
mergedData <- net2KMEs %>% tibble::rownames_to_column(var = "CGid") %>% dplyr::select(CGid) %>% mutate(net2 = net2mergedColors) %>% left_join(mergedData) %>% mutate(net2Matched = as.vector(matchLabels(net2, reference = net1))) %>% filter(CGid %in% rownames(net2KMEs)) 
identical(rownames(net2KMEs),mergedData$CGid) 

net2mergedColors <- mergedData$net2Matched

newModulName <- mergedData %>% dplyr::select(net2, net2Matched) %>% distinct() %>% dplyr::rename(modules =net2, net1Modules = net2Matched)


colnames(net2MEs) <- sapply(colnames(net2MEs), function(x){
  a  <- newModulName$net1Modules[grepl(paste("^",gsub("ME", "", x), "$", sep = ""), newModulName$modules)]
  a  <- paste("ME", a, sep = "")
})

colnames(net2KMEs) <- sapply(colnames(net2KMEs), function(x){
  a  <- newModulName$net1Modules[grepl(paste("^",gsub("MM.", "", x), "$", sep = ""), newModulName$modules)]
  a  <- paste("MM.", a, sep = "")
})

colnames(net2KMEs)[grepl("MM.$", colnames(net2KMEs))] <- "modules"

net2KMEs <- net2KMEs %>% tibble::rownames_to_column(var = "CGid") %>% left_join(mergedData, by = "CGid") %>% dplyr::select(-net2, -net1, -modules)%>% mutate(net2Matched = as.vector(net2Matched)) %>% dplyr::rename(modules = net2Matched)%>% tibble::column_to_rownames(var = "CGid")


# lapply(list("net2MEs", "net2KMEs", "net2mergedColors"), function(x){
#   a <- get(paste(x))
#   saveRDS(a, paste("WGCNA results/All, mappability filter/", x, "_matched", ".RDS", sep = ""))
# })

```


```{r cons networks}
#
load("~/Google Drive/My Drive/Amin documents/Steve projects/Research projects/WGCNA of the total data/WGCNA results/Consensus WGCNA 57 species tissues/Consensus-NetworkConstruction-mammals.RData")
cNetSpeciesTissuesColors <- moduleColors
cNetSpeciesTissuesTree <- consTree
rm(moduleColors, moduleLabels, consTree, consMEs)

#
load("~/Google Drive/My Drive/Amin documents/Steve projects/Research projects/WGCNA of the total data/WGCNA results/Consensus WGCNA 35 species/Consensus-NetworkConstruction-mammals.RData")
cNetSpeciesColors <- moduleColors
cNetSpeciesTree <- consTree
rm(moduleColors, moduleLabels, consTree, consMEs)

#
load("~/Google Drive/My Drive/Amin documents/Steve projects/Research projects/WGCNA of the total data/WGCNA results/Consensus WGCNA 15 tissues/Consensus-NetworkConstruction-mammals.RData")
cNetTissuesColors <- moduleColors
cNetTissuesTree <- consTree
rm(moduleColors, moduleLabels, consTree, consMEs)

# 

colorFiles <- list.files("WGCNA results/Tissue specific Consensus network/", pattern = "moduleColors.RDS", full.names = T, recursive = T)
treeFiles <- list.files("WGCNA results/Tissue specific Consensus network/", pattern = "consTree.RDS", full.names = T, recursive = T)
targets <- gsub("(_moduleColors.RDS)|(AllRegions)", "", list.files("WGCNA results/Tissue specific Consensus network/", pattern = "moduleColors.RDS", full.names = F, recursive = F))


tissueConsColors <- bind_cols(lapply(1:length(colorFiles), function(x){
  a <- readRDS(colorFiles[x])
}))
names(tissueConsColors) <- paste("cNetSpecies",targets, "Colors", sep = "")

tissueConsTrees <- lapply(1:length(treeFiles), function(x){
  a <- readRDS(treeFiles[x])
})
names(tissueConsTrees) <- paste("cNetSpecies",targets, "Tree", sep = "")

```

```{r combined matched colors}
matchedNetworks <- data.frame(CGid = rownames(net1KMEs),net1 = net1mergedColors, 
                    cNet57SpeciesTissues = as.vector(matchLabels(cNetSpeciesTissuesColors, reference = net1mergedColors)),
                    cNet35Species = as.vector(matchLabels(cNetSpeciesColors, reference = net1mergedColors)),
                    cNet15Tissues = as.vector(matchLabels(cNetTissuesColors, reference = net1mergedColors)), 
                    cNet27SpeciesBlood = as.vector(matchLabels(tissueConsColors$cNetSpeciesBloodColors, 
                                                               reference = net1mergedColors)),
                    cNet7SpeciesBrain= as.vector(matchLabels(tissueConsColors$cNetSpeciesBrainColors, 
                                                               reference = net1mergedColors)),
                    cNet10SpeciesLiver = as.vector(matchLabels(tissueConsColors$cNetSpeciesLiverColors, 
                                                               reference = net1mergedColors)),
                    cNet30SpeciesSkin = as.vector(matchLabels(tissueConsColors$cNetSpeciesSkinColors, 
                                                               reference = net1mergedColors))) %>% left_join(mergedData[,c("CGid", "net2Matched")]) %>% dplyr::rename(net2 = net2Matched)

#saveRDS(matchedNetworks, "all networks with matched colors.RDS")
matchedNetworks <- readRDS("all networks with matched colors.RDS")
```


```{r EWAS results}
EWAS_Age <- read.delim("~/Dropbox/HorvathLabCoreMembers/Ake/UniversalClock/EWAS/final/EWASresult/Metal_pgm5_combine_all_species_tissue_stouffer_step2_1.HG38.txt.gz") %>% dplyr::select("CpG", "Meta.Z") %>% dplyr::rename(CGid = CpG, EWAS_Age = Meta.Z)

EWAS_Age_Blood <- read.csv("~/Dropbox/HorvathLabCoreMembers/Ake/UniversalClock/EWAS/final/EWASresult/Metal_pgm6_combine_all_species_Blood_tissue_stouffer_1.HG38.csv.gz") %>% dplyr::select("CpG", "Meta.Z") %>% dplyr::rename(CGid = CpG, EWAS_Age_blood = Meta.Z)

files <- list.files("~/Dropbox/HorvathLabCoreMembers/Ake/UniversalClock/EWAS/final/EWASresult/")
files <- files[grep(".gz$", files)]
targetFiles <- c("all_species_tissue", "Blood","Brain","Cortex", "Liver", "Skin", "Muscle")
targets <- c("all", "blood","brain","cortex", "liver", "skin", "muscle")

all <- read.delim("~/Dropbox/HorvathLabCoreMembers/Ake/UniversalClock/EWAS/final/EWASresult/Metal_pgm5_combine_all_species_tissue_stouffer_step2_1.HG38.txt.gz")


targetList <- lapply(2:length(targetFiles), function(x){
  i <- grep(targetFiles[x], files)
  r <- read.csv(file = paste("~/Dropbox/HorvathLabCoreMembers/Ake/UniversalClock/EWAS/final/EWASresult/", files[i], sep = ""))
})

names(targetList) <- targets[2:7]
targetList$all <- all
targetList <- targetList[c(7,1:6)]

EWAS_Age <- plyr::join_all(lapply(1:length(targetList), function(x){
  nam <- paste("EWAS_age", names(targetList)[x], sep = "_")
  a <- targetList[[x]] %>% dplyr::select("CpG", "Meta.Z") %>% setnames(new= c("CGid", paste(nam)))
}))

```

```{r EWAS max age}
### Use Pvalue and Correlations to get Z values:
## NOTE: change mypath accordingly to your path, and if it's Eutherian path, change them to Eutherian path
generic_association = readRDS("~/Steve Horvath Lab Dropbox/Amin Haghani/HorvathLabCoreMembers/Caesar/ProjectEWASmaxlifespan/allEWAS_maxlifespan/Feb2022/Eutherians_allEWAS_correlation.RDS")# corrs file
myanalysis = "lifespan" ## if this is phylo results, please change to "phylo", 
# in the next step, getStouffer function will summarise the columns to those column tissues accordingly

### convert correlation to Z
stacked.p = readRDS("~/Steve Horvath Lab Dropbox/Amin Haghani/HorvathLabCoreMembers/Caesar/ProjectEWASmaxlifespan/allEWAS_maxlifespan/Feb2022/Eutherians_allEWAS_pvalue.RDS") # pvalue file
generic_association = sign(generic_association) * (- qnorm(stacked.p/2))


# If you want to summarise other columns, e.g. no logmaxAgeCaesar, but just maxAgeCaesar, or only a few tissues, 
# you can change mycols = c(...) accordingly
# Gneeric Tissue columns, (in order): All, blood, skin, liver, brain, muscle.  Phylo Tissues: All, Blood, skin, Liver, Brain, Muscle 
if(myanalysis == "generic") {
  mycols = c(1, 8, 15, 22, 25, 28)
} else {
  mycols = c(1, 4, 7, 12, 15, 18)
}

getStouffer <- function(mydf, myweights = NA) {
  if(!length(myweights) == ncol(mydf)) {
    stop("Number of weights do not match number of columns")
  }
  return(apply(mydf, 1, function(x) return(sum(x * myweights) / sqrt(sum(myweights^2)))))
}
ns = sapply(strsplit(colnames(generic_association)[mycols], "\\.N"), "[", c(2))
ns = as.numeric(ns)[-1]
tissuesStouffer = getStouffer(generic_association[, mycols[-1]], myweights = sqrt(ns))
#stacked = cbind(stacked, tissuesStouffer); colnames(stacked)[ncol(stacked)] = "tissuesStouffer"
generic_association <- as.data.frame(generic_association)
generic_association$Meta <- tissuesStouffer
## Code produces a new column named "tissuesStouffer"
####

Log <- grep("(Log)", colnames(generic_association))
OrderALL <- grep("(OrderALL)|(OrderAll)", colnames(generic_association))
Log <- Log[Log%in%OrderALL]

targets <- c("ALL","Blood", "Brain","Liver","Skin","Meta")
i <- sapply(targets, function(x){
  if(x=="Meta"){
   a <-  grep("Meta", colnames(generic_association))
  } else {
    a <- Log[which(Log %in% grep(paste("Tissue", x, sep = ""), colnames(generic_association)))]
  }
  return(a)
})

EWAS_maxAge <- generic_association[,i] %>% setnames(new = paste("EWAS_lifespan", tolower(targets), sep = "_")) %>% tibble::rownames_to_column(var = "CGid")

rm(generic_association, stacked.p)
```

```{r phylo}
### Use Pvalue and Correlations to get Z values:
## NOTE: change mypath accordingly to your path, and if it's Eutherian path, change them to Eutherian path
phylogenetic_association <- readRDS("~/Steve Horvath Lab Dropbox/Amin Haghani/HorvathLabCoreMembers/Caesar/ProjectEWASmaxlifespan/allEWAS_maxlifespan/Feb2022/Eutherians_allPhyloEWAS_zvalue.RDS")

# corrs file
myanalysis = "phylo" ## if this is phylo results, please change to "phylo", 
# in the next step, getStouffer function will summarise the columns to those column tissues accordingly


# If you want to summarise other columns, e.g. no logmaxAgeCaesar, but just maxAgeCaesar, or only a few tissues, 
# you can change mycols = c(...) accordingly
# Gneeric Tissue columns, (in order): All, blood, skin, liver, brain, muscle.  Phylo Tissues: All, Blood, skin, Liver, Brain, Muscle 
if(myanalysis == "generic") {
  mycols = c(1, 8, 15, 22, 25, 28)
} else {
  mycols = c(1, 4, 7, 12, 15, 18)
}

getStouffer <- function(mydf, myweights = NA) {
  if(!length(myweights) == ncol(mydf)) {
    stop("Number of weights do not match number of columns")
  }
  return(apply(mydf, 1, function(x) return(sum(x * myweights) / sqrt(sum(myweights^2)))))
}
ns = sapply(strsplit(colnames(phylogenetic_association)[mycols], "\\.N"), "[", c(2))
ns = as.numeric(ns)[-1]
tissuesStouffer = getStouffer(phylogenetic_association[, mycols[-1]], myweights = sqrt(ns))
#stacked = cbind(stacked, tissuesStouffer); colnames(stacked)[ncol(stacked)] = "tissuesStouffer"
phylogenetic_association <- as.data.frame(phylogenetic_association)
phylogenetic_association$Meta <- tissuesStouffer

## Code produces a new column named "tissuesStouffer"
####


Log <- grep("(Log)", colnames(phylogenetic_association))
OrderALL <- grep("(OrderALL)|(OrderAll)", colnames(phylogenetic_association))
Log <- Log[Log%in%OrderALL]

i <- sapply(targets, function(x){
  if(x=="Meta"){
   a <-  grep("Meta", colnames(phylogenetic_association))
  } else {
    a <- Log[which(Log %in% grep(paste("Tissue", x, sep = ""), colnames(phylogenetic_association)))]
  }
  return(a)
})


EWAS_maxAgePhylo <- phylogenetic_association[,i] %>% setnames(new = paste("EWAS_lifespanPhylo", tolower(targets), sep = "_")) %>% tibble::rownames_to_column(var = "CGid")


rm(phylogenetic_association, stacked.p)

```


```{r traits}
# bValsNoMars <- readRDS("DNAmDataNoMarsupials.RDS")

# samplesNoMars <- readRDS("samplesNoMars.RDS")

datExpr <- t(bValsNoMars)
##
orders <- samplesNoMars %>% group_by(Order) %>% tally() %>% filter(n>100)
orders <- orders$Order

Order <- as.data.frame(sapply(orders, function(x){
  a <- samplesNoMars %>% mutate(x = ifelse(Order==paste(x), 1, 0)) %>% dplyr::select(x)
  return(a)
}))
colnames(Order) <- orders


datTraits <-  samplesNoMars  %>% dplyr::select(Basename, Female,relativeAge, maximum_age, average_weight, Gestation.days, Age.SexualMaturity)%>% mutate(Female = factor(Female, levels = c(0,1))) %>% droplevels() %>% mutate(log_max_age = log(maximum_age), log_ave_weight = log(average_weight)) %>% dplyr::select(Basename, Female, relativeAge,log_max_age,log_ave_weight, Gestation.days, Age.SexualMaturity) 

datTraits <- samplesNoMars %>% filter(!is.na(relativeAge)&!is.na(Female))%>% mutate(log_max_age = log(maximum_age), log_ave_weight = log(average_weight)) %>% mutate(adjustedRelativeAge = residuals(lm(relativeAge~Tissue+SpeciesLatinName+Female))) %>% mutate(adjustedMaxAge = residuals(lm(log_max_age~relativeAge+Tissue+Female)))%>% mutate(adjustedWeight = residuals(lm(log_ave_weight~relativeAge+Tissue+Female))) %>% dplyr::select(Basename, adjustedRelativeAge, adjustedMaxAge, adjustedWeight) %>% right_join(datTraits)  %>% tibble::column_to_rownames(var = "Basename") %>% relocate(adjustedRelativeAge, .after=relativeAge) %>% relocate(adjustedMaxAge, .after=log_max_age)%>% relocate(adjustedWeight, .after=log_ave_weight)

datTraits <- datTraits[samplesNoMars$Basename,]%>% bind_cols(Order) 


traits <- corAndPvalue(datExpr,datTraits)$cor

# traits <- as.data.frame(traits) %>% tibble::rownames_to_column(var = "CGid") %>% left_join(EWAS_Age) %>% left_join(EWAS_maxAge) %>% left_join(EWAS_maxAgePhylo)  %>%tibble::column_to_rownames(var = "CGid") %>% relocate(starts_with("EWAS_age"), .after=adjustedRelativeAge) %>% relocate(starts_with("EWAS_lifespan"), .after=Age.SexualMaturity)

traits <- as.data.frame(traits) %>% tibble::rownames_to_column(var = "CGid")  %>% left_join(EWAS_maxAge) %>% left_join(EWAS_maxAgePhylo)  %>%tibble::column_to_rownames(var = "CGid") %>% relocate(starts_with("EWAS_lifespan"), .after=Age.SexualMaturity)

names(traits) <- gsub("(max_age)|(maximum_age)|(MaxAge)", "lifespan",names(traits))

traits <- traits %>% dplyr::select(-ends_with("meta"))
traits.Color=numbers2colors(traits,signed=T)
names(traits.Color) <- names(traits)
```

```{r dendogram figures}
sum2 <- readxl::read_xlsx("Drafts/revision 2/7- Supplementary data, Haghani et al, V9.xlsx", sheet = 3) %>% mutate(moduleColors=gsub("ME", "",Modules )) %>% mutate(group=Group)


colnames(matchedNetworks)[2:12] <-c("Net1 (Eutherians)", "Net2 (Mammalian)", "netEPIC", "net450k", "cNet3.omit.Species+Tissues", 
                                 "cNet4.omit.Species", "cNet5.omit.Tissues", "cNet6.Blood.omit.Species","cNet7.Brain.omit.Species", "cNet8.Liver.omit.Species", "cNet9.Skin.omit.Species") 

i <- sapply(c("Net1 \\(Eutherians\\)", "Net2 \\(Mammalian\\)", "cNet3.omit.Species\\+Tissues", 
                                 "cNet4.omit.Species", "cNet5.omit.Tissues", "cNet6.Blood.omit.Species","cNet7.Brain.omit.Species", "cNet8.Liver.omit.Species", "cNet9.Skin.omit.Species"), function(x){grep(x, colnames(matchedNetworks))})

j <- sapply(c("relativeAge", "log_lifespan", "log_ave_weight", "EWAS_lifespanPhylo_blood", "EWAS_lifespanPhylo_brain", "EWAS_lifespanPhylo_liver", "EWAS_lifespanPhylo_skin", "2.Proboscidea"), function(x){grep(x, colnames(traits))})

labels <- sum2 %>% dplyr::select(moduleColors, group)%>% mutate(labels=1:55) %>% mutate(labels = ifelse(group=="unclear", NA, labels)) %>% right_join(dplyr::select(.data=matchedNetworks, "Net1 (Eutherians)", CGid), by=c("moduleColors"="Net1 (Eutherians)")) %>% tibble::column_to_rownames("CGid")%>% mutate(moduleColors = ifelse(group=="unclear", "", moduleColors))
labels <- labels[matchedNetworks$CGid,]

dat <- cbind(matchedNetworks[,i], traits.Color[,-j])

pdf(file = "dendgram Network1.pdf", wi = 5, he = 12)
plotDendroAndColors(net1geneTree, cbind(matchedNetworks[,i], traits.Color[,-j]),
                    c(colnames(matchedNetworks)[i], colnames(traits)[-j]),
                    dendroLabels = F, hang = 0.03,
                    addGuide = T,guideHang = 0.05, main="Network 1",marAll = c(1,10,1,1), cex.colorLabels = 0.8, colorHeightMax = 0.8, cex.rowText = 0.7, rowText = spaste(labels$moduleColors), addTextGuide = T, textPositions = 1, rowWidths = c(1, 8, rep(0.8, 27)))
dev.off()

allTrees <- append(list(net1geneTree,net2geneTree,netEPICgeneTree, net450geneTree, cNetSpeciesTissuesTree, cNetSpeciesTree, cNetTissuesTree), tissueConsTrees)


pdf(file = paste("dendogram", "all", ".pdf"), width = 5, height = 5)
par(mfrow=c(3,3))
for(x in c(2,6:12)){
plotDendroAndColors(allTrees[[x-1]], cbind(matchedNetworks[,c(2,3,6:12)], traits.Color[,-j]),
                    c(colnames(matchedNetworks)[c(2,3,6:12)], colnames(traits)[-j]),
                    dendroLabels = FALSE, hang = 0.001,
                    addGuide = F, guideHang = 0.001, main=paste(names(matchedNetworks)[x]),marAll = c(1,5,1,1), cex.colorLabels = 0.5, colorHeightMax = 0.8)
}
dev.off()

```