Sickle Cell Disease - Weighted Gene Co-Expression Network Analysis (WGCNA)
Yann Ilboudo, Melanie Garrett, Allison Ashley-Koch, Marilyn Telen, Guillaume Lettre
April 25th 2018
1 Abstract
This markdown document is meant to illustrate all the results presented in Ilboudo., et al (2018), all the preprocessing done on the data and the generation of tables as well as figures. The document is interactive and can be be easily navigated moving from one section to another.
Background: Sickle cell disease (SCD) is a monogenic disease caused by mutations in the β-globin gene. The complications related to the disease are systemic as they impact multiple organ systems. Our goal in this study was to identify metabolome changes contributing to SCD-related severity.
Employing both targeted and untargeted approaches, we profiled the plasma of 706 SCD patients using liquid chromatography tandem mass spectrometry. The cohort included 406 French patients of recent African descent and 300 African Americans from the southeastern US. In total, we measured the levels of 61 known and 2,100 unknown metabolites. We applied weighted gene correlation network analysis (WGCNA) algorithms to account for correlations among metabolites and identify specific metabolomic modules associated with SCD-related complications
Several chuncks of code for this analysis are adapted from WGCNA tutorial.
2 Set up
Load or install libraries necessary to perform analysis
source("http://bioconductor.org/biocLite.R")
biocLite(c("AnnotationDbi", "impute", "GO.db", "preprocessCore","WGCNA"))
install.packages(pkgs=c("WGCNA", "flashClust", "ggplot2",
"limma","factoextra","devtools","Biobase","sva","dplyr","magrittr"))
devtools::install_github("yilboudo/myUtils")setwd('~/Desktop/Desktop - MacBook Pro/Lab_Stuff/Metabolite_All_Data_dec52017/wgcna_scd/data/')2.1 Load the libraries and custom scripts
rm(list=ls())
options(stringsAsFactors = F)
library(myUtils)
library(limma)
library(WGCNA)
library(factoextra)
library(devtools)
library(Biobase)
library(sva)
library(dplyr)
library(gdata)
library(data.table)2.2 Load actual datasets
allsamples_new <- read.table("~/Desktop/Desktop - MacBook Pro/Lab_Stuff/Metabolite_All_Data_dec52017/wgcna_scd/data/GENMOD_OMG_metabolites_no_batch_correction.txt",sep="\t",h=T)
batches <- read.table("~/Desktop/Desktop - MacBook Pro/Lab_Stuff/Metabolite_All_Data_dec52017/wgcna_scd/data/GENMOD_OMG_batch_information.txt",sep="\t",h=T)2.3 Align datasets such that metabolites file and batch information are in the same order
#Load dplyr
slice <- dplyr::slice
#align data in the right order
x <- rownames(allsamples_new)
batches <- batches %>% slice(match(x, CLIENT_IDS))
batches <- as.data.frame(batches)2.4 Display cohorts by batches
#image_name1 = 'all_cohorts_batches_pre_correction.pdf'
#pdf(file = image_name1, width = 10, height = 10)
res.pca_uncorrected <- prcomp(allsamples_new, scale = T)
fviz_pca_ind(res.pca_uncorrected, label="none", habillage=batches$Batches,
addEllipses=F, ellipse.level=0.95, palette = "igv") #+ xlim(-50,70) + ylim(-50,50)
#dev.off()
2.5 Apply batch effect correction with ComBat
allsamples_new_combat <- run_combat(allsamples_new, batches, 2)
res.pca_corrected <- prcomp(allsamples_new_combat, scale = TRUE)
#save resulting output
allsamples_new_combat$CLIENT_IDs <- batches$CLIENT_IDS
write.table(allsamples_new_combat, file = "GENMOD_OMG_metabolites_with_batch_correction.txt",
quote = FALSE, row.names = FALSE, col.names = TRUE,sep="\t")2.6 Display cohorts after batch effect correction
image_name2 = 'all_cohorts_batches_post_correction.pdf'
pdf(file = image_name2, width = 10, height = 10)
fviz_pca_ind(res.pca_corrected, label = "none", habillage=batches$Batches,
addEllipses=F, ellipse.level=0.95, palette = "igv") + xlim(-65,40) + ylim(-65,65)
dev.off()
3 Running weighted co-expression analysis
#Load dataset
#Metabolites
allsamples_new_combat <- read.table("~/Desktop/Desktop - MacBook Pro/Lab_Stuff/Metabolite_All_Data_dec52017/wgcna_scd/data/GENMOD_OMG_metabolites_with_batch_correction.txt",sep="\t",h=T,stringsAsFactors = F)
#Phenotypes
large_data_genmod_duke_n705 <-read.table("~/Desktop/Desktop - MacBook Pro/Lab_Stuff/Metabolite_All_Data_dec52017/wgcna_scd/data/genmod_duke_phenotypes_n705_v2.txt",sep="\t",h=T,stringsAsFactors = F)
#Clean metabolite data
rownames(allsamples_new_combat) <- allsamples_new_combat$CLIENT_IDs
allsamples_new_combat <- allsamples_new_combat[,-c(which(names(allsamples_new_combat)=="CLIENT_IDs"))]
#Apply Cyclic Loess Normalization
allsamples_new_combat_norm <- normalizeCyclicLoess(allsamples_new_combat, weights = NULL, span=0.7, iterations = 3, method = "fast")
allsamples_new_combat_norm_cl <- as.data.frame(allsamples_new_combat_norm)
#Apply clustering
sampleTree = hclust(dist(allsamples_new_combat_norm_cl), method = "average")3.1 Step 1. Display sample dendogram - looking for outliers
sizeGrWindow(12,9)
#pdf(file = "allsamples_Clustering.pdf", width = 12, height = 9)
#par(cex = 0.6);
#par(mar = c(0,4,2,0))
plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5,
cex.axis = 1.5, cex.main = 2)
abline(h = 52, col = "red")
#dev.off()
Re-cluster samples and remove outliers from dataset
# Determine cluster under the line
clust = cutreeStatic(sampleTree, cutHeight = 52, minSize = 10)
keepSamples = (clust==1)
allsamples_norm_cl_trimmed = allsamples_new_combat_norm_cl[keepSamples, ]
nMetabo = ncol(allsamples_norm_cl_trimmed)
nSamples = nrow(allsamples_norm_cl_trimmed)Remove outliers from phenotype dataset
allTraits = large_data_genmod_duke_n705;
datTraits <- allTraits[match(rownames(allsamples_norm_cl_trimmed), allTraits$CLIENT_ID),]
#dim(datTraits)
collectGarbage();Cluster samples showing phenotype information on blood traits and complications
sampleTree2 = hclust(dist(allsamples_norm_cl_trimmed), method = "average")
#Convert traits to a color representation: white means low, red means high, grey means missing entry
blood_traitColors = numbers2colors(datTraits[,c(which(names(datTraits)=="eGFR"),which(names(datTraits)=="Hematocrit"),
which(names(datTraits)=="HemoglobinF"),which(names(datTraits)=="MCH"),
which(names(datTraits)=="MCV"),which(names(datTraits)=="Platelets"),
which(names(datTraits)=="RBC"),which(names(datTraits)=="Reticulocyte"),
which(names(datTraits)=="Whitebloodcount"),which(names(datTraits)=="Bilirubin"),which(names(datTraits)=="MCHC"),
which(names(datTraits)=="LDH"))],signed=T)
#blood_traitColors = suppressWarnings(numbers2colors(datTraits[,c(69,27,24)],signed=T))
# Convert traits to a color representation: white means low, red means high, grey means missing entry
complications_traitColors = numbers2colors(datTraits[,c(which(names(datTraits)=="Legulcer"),which(names(datTraits)=="Cholecystectomy"), which(names(datTraits)=="AscepticNecrosis"),which(names(datTraits)=="Priapism"),
which(names(datTraits)=="Retinopathy"),which(names(datTraits)=="Survival"))],signed=T)Display graphics for blood traits and for complications
#pdf(file = "allsamples_with_bloodtraits_clustering2.pdf", width = 20, height = 15);
plotDendroAndColors(sampleTree2, blood_traitColors,
groupLabels = names(datTraits[,c(which(names(datTraits)=="eGFR"),which(names(datTraits)=="Hematocrit"),
which(names(datTraits)=="HemoglobinF"),which(names(datTraits)=="MCH"),
which(names(datTraits)=="MCV"),which(names(datTraits)=="Platelets"),
which(names(datTraits)=="RBC"),which(names(datTraits)=="Reticulocyte"),
which(names(datTraits)=="Whitebloodcount"),which(names(datTraits)=="Bilirubin"),
which(names(datTraits)=="MCHC"),
which(names(datTraits)=="LDH"))]), main = "Sample dendrogram and trait heatmap")
#dev.off()
#pdf(file = "allsamples_with_complication_traits_clustering.pdf", width = 20, height = 15);
plotDendroAndColors(sampleTree2, complications_traitColors,
groupLabels = names(datTraits[,c(which(names(datTraits)=="Legulcer"),
which(names(datTraits)=="Cholecystectomy"),
which(names(datTraits)=="AscepticNecrosis"),which(names(datTraits)=="Priapism"),
which(names(datTraits)=="Retinopathy"),which(names(datTraits)=="Survival"))]),
main = "Sample dendrogram and trait heatmap")
#dev.off()
Save resulting data frames
save(allsamples_norm_cl_trimmed, datTraits, file = "allsamples-01-dataInput.RData")3.2 Step 2. Determine levels of connectivity at different powers
powers = c(c(1:10), seq(from = 12, to=20, by=2))
sft = pickSoftThreshold(allsamples_norm_cl_trimmed,corFnc = cor, networkType = 'signed',powerVector = powers, verbose = 5)Choosing the soft-thresholding power: analysis of network topology
#pdf(file = "allsamples_power.pdf", width = 12, height = 9)
par(mfrow = c(1,2));
cex1 = 0.8;
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
main = paste("Scale independence"));
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
labels=powers,cex=cex1,col="red");
# this line corresponds to using an R^2 cut-off of h
abline(h=0.80,col="red")
plot(sft$fitIndices[,1], sft$fitIndices[,5],
xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
main = paste("Mean connectivity"))
text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")
#dev.off()
3.3 Step 3. Module Construction and Identification
#Co-expression similarity and adjacency
softPower = 6;
adjacency = adjacency(allsamples_norm_cl_trimmed, power = softPower);
#Topological Overlap Matrix (TOM)
TOM = TOMsimilarity(adjacency);
dissTOM = 1-TOM
#Clustering using TOM
geneTree = hclust(as.dist(dissTOM), method = "average");# Plot the resulting clustering tree (dendrogram)
pdf(file = "allsamples_metabo_networks.pdf", width = 12, height = 9)
plot(geneTree, xlab="", sub="", main = "Metabolite clustering on TOM-based dissimilarity",
labels = FALSE, hang = 0.04);
dev.off()Module identification using dynamic tree cut
minModuleSize = 7;
# Module identification using dynamic tree cut:
dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM,
deepSplit = 2, pamRespectsDendro = FALSE,
minClusterSize = minModuleSize);
table(dynamicMods)
# Convert numeric lables into colors
dynamicColors = labels2colors(dynamicMods)
table(dynamicColors)Plot Dendogram with colors underneath
# Plot the dendrogram and colors underneath
sizeGrWindow(8,6)
pdf(file = "allsamples_metabo_networks_color.pdf", width = 12, height = 9)
plotDendroAndColors(geneTree, dynamicColors, "Dynamic Tree Cut",
dendroLabels = FALSE, hang = 0.03,
addGuide = TRUE, guideHang = 0.05,
main = "Metabolite dendrogram and module colors")
dev.off()
##Step 4. Calculate the module eigengenes and combine closely related modules
# Calculate eigengenes
MEList = moduleEigengenes(allsamples_norm_cl_trimmed, colors = dynamicColors)
MEs = MEList$eigengenes
# Calculate dissimilarity of module eigengenes
MEDiss = 1-cor(MEs);
# Cluster module eigengenes
METree = hclust(as.dist(MEDiss), method = "average");
# Select threshold for merging modules together
MEDissThres = 0.25
# Plot the cut line into the dendrogram
#abline(h=MEDissThres, col = "red")
# Call an automatic merging function
merge = mergeCloseModules(allsamples_norm_cl_trimmed, dynamicColors, cutHeight = MEDissThres, verbose = 3)
# The merged module colors
mergedColors = merge$colors;
# Eigengenes of the new merged modules:
mergedMEs = merge$newMEs;Display new modules
#pdf(file = "allsamples_metabo_dendo.pdf", width = 12, height = 9)
plotDendroAndColors(geneTree, cbind(dynamicColors, mergedColors),
c("Dynamic Tree Cut", "Merged dynamic"),
dendroLabels = FALSE, hang = 0.03,
addGuide = TRUE, guideHang = 0.05)
#dev.off()
Rename module colors and save results
# Rename to moduleColors
moduleColors = mergedColors
# Construct numerical labels corresponding to the colors
colorOrder = c("grey", standardColors(50));
moduleLabels = match(moduleColors, colorOrder)-1;
MEs = mergedMEs;
# Save module colors and labels for use in subsequent parts
save(MEs, moduleLabels, moduleColors, geneTree, file = "allsamples-02-networkConstruction-stepByStep.RData")3.4 Step 5.Remove modules with poor preservation scores, and re-draw topological overlap map (TOM) and eigengenes
lnames3 = load(file = "~/Desktop/Desktop - MacBook Pro/Lab_Stuff/Metabolite_All_Data_dec52017/module_preservation_stats.RData")
preserved_modules <- preservations_stats_all_df[preservations_stats_all_df$V2 > 10,]
restGenes= (moduleColors %in% preserved_modules$V1) & !(moduleColors %in% "grey") & !(moduleColors %in% "gold")
SubGeneNames <- colnames(allsamples_norm_cl_trimmed)
dissTOM = 1-TOMsimilarityFromExpr(allsamples_norm_cl_trimmed[,restGenes], power = 6);
# Transform dissTOM with a power to make moderately strong connections more visible in the heatmap
plotTOM = dissTOM^10;
colnames(dissTOM) =rownames(dissTOM) =SubGeneNames[restGenes]
hier1=hclust(as.dist(dissTOM), method="average" )
# Set diagonal to NA for a nicer plot
diag(plotTOM) = NA;
# Call the plot function
#sizeGrWindow(9,9)
#pdf(file = "all_samples_network_plot2.pdf", width = 10, height = 10)
TOMplot(plotTOM, hier1, moduleColors[restGenes], main = "Network heatmap plot, all genes")
#dev.off()
MEList = moduleEigengenes(allsamples_norm_cl_trimmed[,restGenes], colors = moduleColors[restGenes])
MEs = MEList$eigengenes
MET=orderMEs(MEs)
#pdf(file = "all_samples_eigengenes_plot2.pdf", width = 10, height = 10)
plotEigengeneNetworks(MET, "", marDendro = c(0,4,1,2), marHeatmap = c(3,4,1,2))
#dev.off()
4 Relate module to external traits
4.1 Relate eigen genes to blood traits and complications
lnames1 = load(file = "~/Desktop/Desktop - MacBook Pro/Lab_Stuff/Metabolite_All_Data_dec52017/allsamples-01-dataInput.RData");
#The variable lnames contains the names of loaded variables.
lnames12 = load(file = "~/Desktop/Desktop - MacBook Pro/Lab_Stuff/Metabolite_All_Data_dec52017/allsamples-03-associations_annotation_external_traits.RData")
lnames13 = load(file = "~/Desktop/Desktop - MacBook Pro/Lab_Stuff/Metabolite_All_Data_dec52017/allsamples-03-associations_external_traits.RData")
# Load network data saved in the second part.
lnames2 = load(file = "~/Desktop/Desktop - MacBook Pro/Lab_Stuff/Metabolite_All_Data_dec52017/allsamples-02-networkConstruction-stepByStep.RData");
datTraits_for_assoc <- datTraits
# Define numbers of metabolites and samples
nGenes = ncol(allsamples_norm_cl_trimmed);
nSamples = nrow(allsamples_norm_cl_trimmed);
# Recalculate MEs with color labels
MEs0 = moduleEigengenes(allsamples_norm_cl_trimmed, moduleColors)$eigengenes
MEs = orderMEs(MEs0)
rownames(MEs) <- rownames(allsamples_norm_cl_trimmed)
#inverse normal transform eigengenes
MEs_for_assoc <- as.data.frame(apply(MEs,2,invnorm))
#Define variable for blood traits, and covariates
blood_traits1 = "eGFR"
cov_names1="SCD_geno,HU"
blood_traits = c("Hematocrit","HemoglobinF","MCH","MCV","MCHC","Platelets","RBC","Reticulocyte","Whitebloodcount", "Bilirubin", "LDH")
cov_names="SCD_geno,HU,Age,Sex"
#Create matrices to store summary statistics values
blood_traits_associations_pval <- matrix(NA,dim(MEs)[2],length(blood_traits))
blood_traits_associations_effect <- matrix(NA,dim(MEs)[2],length(blood_traits))
blood_traits_associations_n <- matrix(NA,dim(MEs)[2],length(blood_traits))
blood_traits_associations_pval1 <- matrix(NA,dim(MEs)[2],length(blood_traits1))
blood_traits_associations_effect1 <- matrix(NA,dim(MEs)[2],length(blood_traits1))
blood_traits_associations_n1 <- matrix(NA,dim(MEs)[2],length(blood_traits1))
rownames(datTraits_for_assoc) <- (datTraits_for_assoc$CLIENT_ID)
#Run association for most traits except eGFR because covariates are not the same
for (i in 1:length(blood_traits)) {
datTraits_trait <- datTraits_for_assoc[complete.cases(datTraits_for_assoc[,which(colnames(datTraits_for_assoc)==blood_traits[i])]),]
keeping_indiv <- rownames(datTraits_trait)
datTraits_trait$Zpheno <- invnorm(datTraits_trait[,which(names(datTraits_trait)==blood_traits[i])])
MEsTrait <- MEs_for_assoc[(rownames(MEs_for_assoc) %in% keeping_indiv),]
covariates <- strsplit(cov_names, split=",")[[1]]
cov=datTraits_trait[,covariates]
observed_pval <- sapply(MEsTrait, linearRegcov_p, y=datTraits_trait$Zpheno, data=cov)
std_err <- sapply(MEsTrait, linearRegcov_stderr, y=datTraits_trait$Zpheno, data=cov)
effect_size <- sapply(MEsTrait, linearRegcov_effect, y=datTraits_trait$Zpheno, data=cov)
sample_size <- sapply(MEsTrait, linearRegcov_samplesize, y=datTraits_trait$Zpheno, data=cov)
blood_traits_associations_pval[,i] <- observed_pval
blood_traits_associations_effect[,i] <- effect_size
blood_traits_associations_n[,i] <- sample_size
}
#Run association for eGFR
for (i in 1:length(blood_traits1)) {
datTraits_trait <- datTraits_for_assoc[complete.cases(datTraits_for_assoc[,which(colnames(datTraits_for_assoc)==blood_traits1[i])]),]
keeping_indiv <- rownames(datTraits_trait)
datTraits_trait$Zpheno <- invnorm(datTraits_trait[,which(names(datTraits_trait)==blood_traits1[i])])
MEsTrait <- MEs_for_assoc[(rownames(MEs_for_assoc) %in% keeping_indiv),]
covariates <- strsplit(cov_names1, split=",")[[1]]
cov=datTraits_trait[,covariates]
observed_pval <- sapply(MEsTrait, linearRegcov_p, y=datTraits_trait$Zpheno, data=cov)
std_err <- sapply(MEsTrait, linearRegcov_stderr, y=datTraits_trait$Zpheno, data=cov)
effect_size <- sapply(MEsTrait, linearRegcov_effect, y=datTraits_trait$Zpheno, data=cov)
sample_size <- sapply(MEsTrait, linearRegcov_samplesize, y=datTraits_trait$Zpheno, data=cov)
blood_traits_associations_pval1[,i] <- observed_pval
blood_traits_associations_effect1[,i] <- effect_size
blood_traits_associations_n1[,i] <- sample_size
}
blood_traits_associations_pval <- cbind(blood_traits_associations_pval, blood_traits_associations_pval1)
blood_traits_associations_effect <- cbind(blood_traits_associations_effect, blood_traits_associations_effect1)
blood_traits_associations_n <- cbind(blood_traits_associations_n, blood_traits_associations_n1)
blood_traits_associations <- as.data.frame(blood_traits_associations_pval)
names(blood_traits_associations) <- c(blood_traits , "eGFR")
rownames(blood_traits_associations) <- names(MEs)
blood_traits_associations_cor <- as.data.frame(blood_traits_associations_effect)
names(blood_traits_associations_cor) <- c(blood_traits , "eGFR")
rownames(blood_traits_associations_cor) <- names(MEs)
#Associations with complications
complications = c("Legulcer","Cholecystectomy","AscepticNecrosis","Priapism", "Retinopathy","Survival")
complications_associations_pval <- matrix(NA,dim(MEs)[2],length(complications))
complications_associations_effect <- matrix(NA,dim(MEs)[2],length(complications))
complications_associations_n <- matrix(NA,dim(MEs)[2],length(complications))
for (i in 1:length(complications)) {
datTraits_trait <- datTraits_for_assoc[complete.cases(datTraits_for_assoc[,which(colnames(datTraits_for_assoc)==complications[i])]),]
keeping_indiv <- rownames(datTraits_trait)
datTraits_trait$Zpheno <- datTraits_trait[,which(names(datTraits_trait)==complications[i])]
MEsTrait <- MEs_for_assoc[(rownames(MEs_for_assoc) %in% keeping_indiv),]
covariates <- strsplit(cov_names, split=",")[[1]]
cov=datTraits_trait[,covariates]
observed_pval <- sapply(MEsTrait, logisticRegcov_p, y=datTraits_trait$Zpheno, data=cov)
effect_size <- sapply(MEsTrait, logisticRegcov_effect, y=datTraits_trait$Zpheno, data=cov)
sample_size <- sapply(MEsTrait, logisticRegcov_samplesize, y=datTraits_trait$Zpheno, data=cov)
complications_associations_pval[,i] <- observed_pval
complications_associations_effect[,i] <- effect_size
complications_associations_n[,i] <- sample_size
}
complications_associations <- as.data.frame(complications_associations_pval)
names(complications_associations) <- complications
rownames(complications_associations) <- names(MEs)
complications_associations_cor <- as.data.frame(complications_associations_effect)
names(complications_associations_cor) <- complications
rownames(complications_associations_cor) <- names(MEs)
pvals <- cbind(complications_associations,blood_traits_associations)
cor <- cbind(complications_associations_cor,blood_traits_associations_cor)
pvals_subset
pvals
pvals_subset <- pvals[c("MEdarkorange",
"MEyellow4",
"MElightcyan",
"MEmediumorchid",
"MEfirebrick4",
"MEmaroon",
"MEdarkgreen",
"MEdarkslateblue",
"MEmagenta",
"MEblue",
"MEbrown4",
"MEskyblue1",
"MEdarkgrey",
"MElightcoral",
"MEsienna3",
"MElightcyan1",
"MEdarkviolet",
"MEdarkorange2",
"MEskyblue",
"MElightsteelblue1",
"MEmediumpurple2",
"MEdarkred",
"MEdarkolivegreen4",
"MEbrown2",
"MEyellowgreen",
"MElightyellow",
"MEbisque4",
"MEroyalblue",
"MElightpink4",
"MEantiquewhite4",
"MEsteelblue",
"MEplum",
"MEgreenyellow",
"MEhoneydew1",
"MEskyblue3",
"MEgrey60",
"MEcyan",
"MEcoral1",
"MEdarkmagenta","MEcoral2","MEbrown","MEplum2","MEgrey"), c("Cholecystectomy", "Hematocrit", "HemoglobinF", "Platelets", "MCH", "MCV", "MCHC", "RBC", "Reticulocyte", "Whitebloodcount", "Bilirubin", "LDH", "eGFR")]
cor_subset <- cor[c("MEdarkorange",
"MEyellow4",
"MElightcyan",
"MEmediumorchid",
"MEfirebrick4",
"MEmaroon",
"MEdarkgreen",
"MEdarkslateblue",
"MEmagenta",
"MEblue",
"MEbrown4",
"MEskyblue1",
"MEdarkgrey",
"MElightcoral",
"MEsienna3",
"MElightcyan1",
"MEdarkviolet",
"MEdarkorange2",
"MEskyblue",
"MElightsteelblue1",
"MEmediumpurple2",
"MEdarkred",
"MEdarkolivegreen4",
"MEbrown2",
"MEyellowgreen",
"MElightyellow",
"MEbisque4",
"MEroyalblue",
"MElightpink4",
"MEantiquewhite4",
"MEsteelblue",
"MEplum",
"MEgreenyellow",
"MEhoneydew1",
"MEskyblue3",
"MEgrey60",
"MEcyan",
"MEcoral1",
"MEdarkmagenta","MEcoral2","MEbrown","MEplum2","MEgrey"), c("Cholecystectomy", "Hematocrit", "HemoglobinF", "Platelets", "MCH", "MCV", "MCHC", "RBC", "Reticulocyte", "Whitebloodcount", "Bilirubin", "LDH", "eGFR")]
MEs_subset <- MEs[,c("MEdarkorange",
"MEyellow4",
"MElightcyan",
"MEmediumorchid",
"MEfirebrick4",
"MEmaroon",
"MEdarkgreen",
"MEdarkslateblue",
"MEmagenta",
"MEblue",
"MEbrown4",
"MEskyblue1",
"MEdarkgrey",
"MElightcoral",
"MEsienna3",
"MElightcyan1",
"MEdarkviolet",
"MEdarkorange2",
"MEskyblue",
"MElightsteelblue1",
"MEmediumpurple2",
"MEdarkred",
"MEdarkolivegreen4",
"MEbrown2",
"MEyellowgreen",
"MElightyellow",
"MEbisque4",
"MEroyalblue",
"MElightpink4",
"MEantiquewhite4",
"MEsteelblue",
"MEplum",
"MEgreenyellow",
"MEhoneydew1",
"MEskyblue3",
"MEgrey60",
"MEcyan",
"MEcoral1",
"MEdarkmagenta","MEcoral2","MEbrown","MEplum2","MEgrey")]Display associations statistics of modules and blood traits and complications in a heatmap
png(file="top_blood_metabo_traits_complications_2022.png", units = "px",width=14250, height=13250, res=350, pointsize=25)
# Will display correlations and their p-values
textMatrix = paste(signif(as.matrix(pvals_subset), 2), "\n(",
signif(as.matrix(cor_subset), 1), ")", sep = "");
dim(textMatrix) = dim(cor_subset)
par(mar = c(7, 11, 2, 3));
# Display the correlation values within a heatmap plot
labeledHeatmap(Matrix = cor_subset[,which(names(pvals_subset)=="Hematocrit"):which(names(pvals_subset)=="eGFR")],
xLabels = names(pvals_subset)[which(names(pvals_subset)=="Hematocrit"):which(names(pvals_subset)=="eGFR")],
yLabels = names(MEs_subset),
ySymbols = names(MEs_subset),
colorLabels = F,
colors = blueWhiteRed(50),
textMatrix = textMatrix[,which(names(pvals_subset)=="Hematocrit"):which(names(pvals_subset)=="eGFR")],
setStdMargins = F,
cex.text = 1.0,
zlim = c(-0.60, 0.40),
main = paste(""))
dev.off()png(file="top_complications_metabo_trait_2022.png", units = "px",width=9000, height=13250, res=350, pointsize=25)
# Will display correlations and their p-values
textMatrix = paste(signif(as.matrix(pvals_subset), 2), "\n(",
signif(as.matrix(cor_subset), 1), ")", sep = "");
dim(textMatrix) = dim(cor_subset)
par(mar = c(7, 10, 2, 3));
# Display the correlation values within a heatmap plot
labeledHeatmap(Matrix = as.matrix(cor_subset[,which(names(pvals_subset)=="Cholecystectomy")]),
xLabels = names(pvals_subset)[which(names(pvals_subset)=="Cholecystectomy")],
yLabels = names(MEs_subset),
ySymbols = names(MEs_subset),
colorLabels = F,
colors = blueWhiteRed(50),
textMatrix = textMatrix[,which(names(pvals_subset)=="Cholecystectomy")],
setStdMargins = F,
cex.text = 1.0,
zlim = c(0, 2),
main = paste(""))
dev.off()Display associations statistics of modules and blood traits and complications in a heatmap
png(file="allsamples_metabo_traits_assoc_blood_final_2022.png", units = "px",width=9000, height=9000, res=350, pointsize=15)
# Will display correlations and their p-values
textMatrix = paste(signif(as.matrix(pvals), 2), "\n(",
signif(as.matrix(cor), 1), ")", sep = "");
dim(textMatrix) = dim(cor)
par(mar = c(7, 8.5, 2, 3));
# Display the correlation values within a heatmap plot
labeledHeatmap(Matrix = cor[,which(names(pvals)=="Hematocrit"):which(names(pvals)=="eGFR")],
xLabels = names(pvals)[which(names(pvals)=="Hematocrit"):which(names(pvals)=="eGFR")],
yLabels = names(MEs),
ySymbols = names(MEs),
colorLabels = F,
colors = blueWhiteRed(50),
textMatrix = textMatrix[,which(names(pvals)=="Hematocrit"):which(names(pvals)=="eGFR")],
setStdMargins = F,
cex.text = 0.7,
zlim = c(-0.60, 0.40),
main = paste("Module-trait relationships"))
dev.off()`
png(file="allsamples_metabo_traits_assoc_complications_final.png", units = "px",width=9000, height=9000, res=350, pointsize=15)
# Will display correlations and their p-values
textMatrix = paste(signif(as.matrix(pvals), 2), "\n(",
signif(as.matrix(cor), 1), ")", sep = "");
dim(textMatrix) = dim(cor)
par(mar = c(7, 8.5, 2, 3));
# Display the correlation values within a heatmap plot
labeledHeatmap(Matrix = cor[,which(names(pvals)=="Legulcer"):which(names(pvals)=="Survival")], #moduleTraitCor[,-c(1,15:19)],
xLabels = names(pvals)[which(names(pvals)=="Legulcer"):which(names(pvals)=="Survival")],
yLabels = names(MEs),
ySymbols = names(MEs),
colorLabels = F,
colors = blueWhiteRed(50),
textMatrix = textMatrix[,which(names(pvals)=="Legulcer"):which(names(pvals)=="Survival")],
setStdMargins = F,
cex.text = 0.7,
zlim = c(0,2),
main = paste("Module-trait relationships"))
dev.off()