WGCNA trait order in heat map
0
1
Entering edit mode
coyoung ▴ 10
@coyoung-17963
Last seen 4.8 years ago

Does anyone know how to edit the order of the traits in the Module−trait relationships heat map or match the order in my clinical file? The order of these traits is way out of wack. Can anyone help?

## Output GlobalNetworkPlots and kMEtable
####################################################################################################################
FileBaseName=paste0(projectFilesOutputTag,"_pwr",power,"_PAMstageTRUE")

pdf(file=paste0(outputfigs,"1.GlobalNetworkPlots-",FileBaseName,".pdf"),width=16,height=12)
#pdf(file=paste0("1.GlobalNetworkPlots-",FileBaseName,".pdf"),width=16,height=12)

## Plot dendrogram with module colors and trait correlations
MEs<-tmpMEs<-data.frame()
MEList = moduleEigengenes(t(cleanDat), colors = net$colors)
MEs = orderMEs(MEList$eigengenes)
colnames(MEs)<-gsub("ME","",colnames(MEs)) #let's be consistent in case prefix was added, remove it.
rownames(MEs)<-rownames(numericMeta)

numericIndices<-unique(c( which(!is.na(apply(numericMeta,2,function(x) sum(as.numeric(x))))), which(!(apply(numericMeta,2,function(x) sum(as.numeric(x),na.rm=T)))==0) ))
#Warnings OK; This determines which traits are numeric and if forced to numeric values, non-NA values do not sum to 0

geneSignificance <- cor(sapply(numericMeta[,numericIndices],as.numeric),t(cleanDat),use="pairwise.complete.obs")
rownames(geneSignificance) <- colnames(numericMeta)[numericIndices]
geneSigColors <- t(numbers2colors(t(geneSignificance),signed=TRUE,lim=c(-1,1),naColor="black"))
rownames(geneSigColors) <- colnames(numericMeta)[numericIndices]

plotDendroAndColors(dendro=net$dendrograms[[1]],
                    colors=t(rbind(net$colors,geneSigColors)),
                    cex.dendroLabels=1.2,addGuide=TRUE,
                    dendroLabels=FALSE,
                    groupLabels=c("Module Colors",colnames(numericMeta)[numericIndices]))

## Plot eigengene dendrogram/heatmap - using bicor
tmpMEs <- MEs #net$MEs
colnames(tmpMEs) <- paste("ME",colnames(MEs),sep="")
MEs[,"grey"] <- NULL
tmpMEs[,"MEgrey"] <- NULL

plotEigengeneNetworks(tmpMEs, "Eigengene Network", marHeatmap = c(3,4,2,2), marDendro = c(0,4,2,0),plotDendrograms = TRUE, xLabelsAngle = 90,heatmapColors=blueWhiteRed(50))


######################
## Find differences between Groups (as defined in Traits input file); Finalize Grouping of Samples for ANOVA




# This gets ANOVA (Kruskal-Wallis) nonparametric p-values for groupwise comparison of interest.
# look at numericMeta (traits data) and choose traits to use for linear model-determination of p value
head(numericMeta)
# Change below line to point to a factored trait, which will define groups for ANOVA
regvars <- data.frame(as.factor( numericMeta$Group ), as.numeric(numericMeta$Age), as.numeric(numericMeta$Sex))
colnames(regvars) <- c("Group","Age","Sex") ## data frame with covaraites incase we want to try multivariate regression
##aov1 <- aov(data.matrix(MEs)~Group,data=regvars) ## ANOVA framework yields same results
lm1 <- lm(data.matrix(MEs)~Group,data=regvars) #sex and age effects are removed by the linear model

pvec <- rep(NA,ncol(MEs))
for (i in 1:ncol(MEs)) {
  f <- summary(lm1)[[i]]$fstatistic ## Get F statistics
  pvec[i] <- pf(f[1],f[2],f[3],lower.tail=F) ## Get the p-value corresponding to the whole model
}
names(pvec) <- colnames(MEs)


######################
## Get sigend kME values
kMEdat <- signedKME(t(cleanDat), tmpMEs, corFnc="bicor")


######################
## Plot eigengene-trait correlations - using p value of bicor for heatmap scale
library(RColorBrewer)
MEcors <- bicorAndPvalue(MEs,numericMeta[,numericIndices])
moduleTraitCor <- MEcors$bicor
moduleTraitPvalue <- MEcors$p


textMatrix = apply(moduleTraitCor,2,function(x) signif(x, 2))
#textMatrix = paste(signif(moduleTraitCor, 2), " (",
#  signif(moduleTraitPvalue, 1), ")", sep = "");
#dim(textMatrix) = dim(moduleTraitCor)
par(mfrow=c(1,1))
par(mar = c(6, 8.5, 3, 3));

## Display the correlation values within a heatmap plot
cexy <- if(nModules>75) { 0.8 } else { 1 }
colvec <- rep("white",1500)
colvec[1:500] <- colorRampPalette(rev(brewer.pal(8,"BuPu")[2:8]))(500)
colvec[501:1000]<-colorRampPalette(c("white",brewer.pal(8,"BuPu")[2]))(3)[2] #interpolated color for 0.05-0.1 p
labeledHeatmap(Matrix = apply(moduleTraitPvalue,2,as.numeric),
               xLabels = colnames(numericMeta)[numericIndices],
               yLabels = paste0("ME",names(MEs)),
               ySymbols = names(MEs),
               colorLabels = FALSE,
               colors = colvec,
               textMatrix = textMatrix,
               setStdMargins = FALSE,
               cex.text = 0.5,
               cex.lab.y= cexy,
               zlim = c(0,0.15),
               main = paste("Module-trait relationships\n bicor r-value shown as text\nHeatmap scale: Student correlation p value"),
               cex.main=0.8)
wgcna WGCNA wgcna package • 1.1k views
ADD COMMENT

Login before adding your answer.

Traffic: 766 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6