Entering edit mode
coyoung
▴
10
@coyoung-17963
Last seen 4.9 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)