title: “Wgcna.trait.module.analysis” author: “Dr UPASNA” date: “2024-06-30” output: html_document
library(flashClust)##
## Attaching package: 'flashClust'## The following object is masked from 'package:stats':
##
## hclustlibrary(WGCNA)## Loading required package: dynamicTreeCut## Loading required package: fastcluster##
## Attaching package: 'fastcluster'## The following object is masked from 'package:flashClust':
##
## hclust## The following object is masked from 'package:stats':
##
## hclust## ##
## Attaching package: 'WGCNA'## The following object is masked from 'package:stats':
##
## corlibrary(WGCNA)
library(tidyr)## Warning: package 'tidyr' was built under R version 4.3.2library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.3.2library(cowplot)## Warning: package 'cowplot' was built under R version 4.3.2library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(stringr)
library(rmarkdown)## Warning: package 'rmarkdown' was built under R version 4.3.3library(WGCNA)
library(GO.db)## Loading required package: AnnotationDbi## Warning: package 'AnnotationDbi' was built under R version 4.3.2## Loading required package: stats4## Loading required package: BiocGenerics## Warning: package 'BiocGenerics' was built under R version 4.3.1##
## Attaching package: 'BiocGenerics'## The following objects are masked from 'package:dplyr':
##
## combine, intersect, setdiff, union## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs## The following objects are masked from 'package:base':
##
## anyDuplicated, aperm, append, as.data.frame, basename, cbind,
## colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
## get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
## match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
## Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
## table, tapply, union, unique, unsplit, which.max, which.min## Loading required package: Biobase## Warning: package 'Biobase' was built under R version 4.3.1## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.## Loading required package: IRanges## Warning: package 'IRanges' was built under R version 4.3.1## Loading required package: S4Vectors## Warning: package 'S4Vectors' was built under R version 4.3.2##
## Attaching package: 'S4Vectors'## The following objects are masked from 'package:dplyr':
##
## first, rename## The following object is masked from 'package:tidyr':
##
## expand## The following object is masked from 'package:utils':
##
## findMatches## The following objects are masked from 'package:base':
##
## expand.grid, I, unname##
## Attaching package: 'IRanges'## The following objects are masked from 'package:dplyr':
##
## collapse, desc, slice##
## Attaching package: 'AnnotationDbi'## The following object is masked from 'package:dplyr':
##
## select###
options(strigsAsFactors = FALSE)
# Uploading data into R and formatting it for WGCNA
cleanDat <- readRDS("~/Desktop/Camk2a.WGCNA.8june/JUNE9.CAMK2A/part2/cleanDat.after.outlier.removal.15samples.rds")
dim(cleanDat)#in this rows is genes and column is samples which we will transpose in WGCNA scripts## [1] 1544 15head(cleanDat, 10)## X.106P2fractionIP X.120P2fractionIP X.117P2fractionIP
## AAK1|Q2M2I8 23.48784 23.92167 23.18945
## AAK1|Q2M2I8.1 25.21936 25.97150 25.12306
## ABAT|P80404 27.55898 27.40838 27.95464
## ABHD12|Q8N2K0 22.68297 23.18231 24.03555
## ABHD16A|O95870 22.20707 19.64799 22.80265
## ABHD6|Q9BV23 19.50613 19.97449 20.11511
## ABI1|Q8IZP0-9 22.68659 20.34861 22.29118
## ABI2|Q9NYB9 22.42038 23.21194 22.08461
## ABLIM1|O14639-6 24.53111 24.49455 20.58515
## ABLIM1|O14639-6.1 22.28356 19.95007 21.25693
## X.123P2fractionIP X.137P2fractionIP X.131P2fractionIP
## AAK1|Q2M2I8 23.31690 23.23492 22.92868
## AAK1|Q2M2I8.1 25.51681 25.18611 25.24855
## ABAT|P80404 27.51978 27.19804 27.35241
## ABHD12|Q8N2K0 23.16306 23.51496 23.37672
## ABHD16A|O95870 22.37072 22.40512 22.94137
## ABHD6|Q9BV23 22.63349 22.86750 19.39300
## ABI1|Q8IZP0-9 22.71380 22.92258 22.31786
## ABI2|Q9NYB9 22.53702 21.89735 22.05724
## ABLIM1|O14639-6 23.79206 24.41973 24.78408
## ABLIM1|O14639-6.1 20.38977 21.76093 21.40035
## X.143P2fractionIP X.133P2fractionIP X.106HomogenateIP
## AAK1|Q2M2I8 23.61881 22.99170 22.93829
## AAK1|Q2M2I8.1 25.07017 25.34287 25.46261
## ABAT|P80404 27.27474 27.16163 27.14351
## ABHD12|Q8N2K0 24.10007 23.79960 23.34190
## ABHD16A|O95870 22.03687 23.16943 22.30628
## ABHD6|Q9BV23 22.65173 23.06062 22.54051
## ABI1|Q8IZP0-9 23.04449 20.06045 22.13664
## ABI2|Q9NYB9 22.25320 19.95309 22.64059
## ABLIM1|O14639-6 20.37294 20.58515 24.66178
## ABLIM1|O14639-6.1 20.24092 21.52648 20.38148
## X.120HomogenateIP X.117HomogenateIP X.123HomogenateIP
## AAK1|Q2M2I8 23.87091 23.85305 23.11254
## AAK1|Q2M2I8.1 25.80004 25.22436 25.56972
## ABAT|P80404 26.87173 26.93971 26.87173
## ABHD12|Q8N2K0 23.36082 22.77758 23.39418
## ABHD16A|O95870 22.79152 22.08915 21.14526
## ABHD6|Q9BV23 19.95007 22.72849 22.30235
## ABI1|Q8IZP0-9 20.65219 22.52323 22.59062
## ABI2|Q9NYB9 22.36293 22.28723 23.15349
## ABLIM1|O14639-6 20.92849 20.70148 19.34104
## ABLIM1|O14639-6.1 20.42097 19.88261 21.65811
## X.131HomogenateIP X.143HomogenateIP X.133HomogenateIP
## AAK1|Q2M2I8 23.40135 19.97449 23.80321
## AAK1|Q2M2I8.1 25.29442 25.52288 25.54158
## ABAT|P80404 27.12597 27.27474 26.92225
## ABHD12|Q8N2K0 22.98165 23.73187 22.74293
## ABHD16A|O95870 22.08461 22.14965 22.18829
## ABHD6|Q9BV23 20.55124 20.20580 21.11443
## ABI1|Q8IZP0-9 22.27947 20.36560 20.19984
## ABI2|Q9NYB9 23.13796 22.92258 22.01720
## ABLIM1|O14639-6 21.25110 20.66485 20.85644
## ABLIM1|O14639-6.1 20.24928 22.51959 19.52090# Manipulate file so it matches the format WGCNA needs here I am transpose the matrix
datExpr <- t(cleanDat)
dim(datExpr) #in this samples are in rows and column as genes## [1] 15 1544#load metadata for trait information
datTraits <- readRDS("~/Desktop/Camk2a.WGCNA.8june/JUNE9.CAMK2A/part2/datTraits.15samples.condition.rds")
datTraits## condition
## X.106P2fractionIP P2fractionIP.CamK2a.LPS
## X.120P2fractionIP P2fractionIP.CamK2a.saline
## X.117P2fractionIP P2fractionIP.CamK2a.LPS
## X.123P2fractionIP P2fractionIP.CamK2a.saline
## X.137P2fractionIP P2fractionIP.CamK2a.LPS
## X.131P2fractionIP P2fractionIP.CamK2a.saline
## X.143P2fractionIP P2fractionIP.CamK2a.LPS
## X.133P2fractionIP P2fractionIP.CamK2a.saline
## X.106HomogenateIP HomogenateIP.CamK2a.LPS
## X.120HomogenateIP HomogenateIP.CamK2a.saline
## X.117HomogenateIP HomogenateIP.CamK2a.LPS
## X.123HomogenateIP HomogenateIP.CamK2a.saline
## X.131HomogenateIP HomogenateIP.CamK2a.saline
## X.143HomogenateIP HomogenateIP.CamK2a.LPS
## X.133HomogenateIP HomogenateIP.CamK2a.salinenumericMeta <- datTraits
table(rownames(datTraits)==rownames(datExpr)) #should return TRUE if datasets align correctly, otherwise your names are out of order##
## TRUE
## 15## Run this to check if there are gene outliers
## 3. Check and Remove Outliers:
sdout=3 #Z.k SD fold for outlier threshold | any samples that are above or below 3 sd will be excluded
outliers.All<-vector()
cleanDat.All<-cleanDat
numericMeta.All<-numericMeta
for (repeated in 1:20) {
normadj <- (0.5+0.5*bicor(cleanDat,use="pairwise.complete.obs")^2)
## Calculate connectivity
netsummary <- fundamentalNetworkConcepts(normadj)
ku <- netsummary$Connectivity
z.ku <- ku-(mean(ku))/sqrt(var(ku))
## Declare as outliers those samples which are more than sdout sd above the mean connectivity based on the chosen measure
outliers <- (z.ku > mean(z.ku)+sdout*sd(z.ku))|(z.ku < mean(z.ku)-sdout*sd(z.ku))
print(paste0("There are ",sum(outliers)," outlier samples based on a bicor distance sample network connectivity standard deviation above ",sdout,". [Round ",repeated,"]"))
targets.All=numericMeta
cleanDat <- cleanDat[,!outliers]
numericMeta <- targets <- targets.All[!outliers,]
outliers.All<-c(outliers.All,outliers)
if (sum(outliers)==0) break;
} #repeat up to 20 times## [1] "There are 0 outlier samples based on a bicor distance sample network connectivity standard deviation above 3. [Round 1]"#All outliers removed #"There are 1 total outlier samples removed in 2 iterations:"
print(paste0("There are ",sum(outliers.All)," total outlier samples removed in ",repeated," iterations:"))## [1] "There are 0 total outlier samples removed in 1 iterations:"#character(0)
names(which(outliers.All)) #:Identified outliar removed!## character(0)# Current pass OL removal output:
#"There are 1 outlier samples based on a bicor distance sample network connectivity standard deviation above 3. [Round 1]"
outliersRemoved<-names(which(outliers.All))
## Enforce <50% missingness (1 less than half of columns (or round down half if odd number of columns))
#Any protein that has less than 50% values within the samples will be dropped
#Leave samples with more than 50% NA's no imputations were done here!
LThalfSamples <- length(colnames(cleanDat)) / 2
LThalfSamples <- LThalfSamples - if ((length(colnames(cleanDat)) %% 2) == 1) {
0.5
} else {
1.0
}
removedRownames <- rownames(cleanDat[which(rowSums(as.matrix(is.na(cleanDat))) > LThalfSamples), ]) # list rows to be removed
if(length(which(rowSums(as.matrix(is.na(cleanDat))) > LThalfSamples))==1) removedRownames<-rownames(cleanDat)[which(rowSums(as.matrix(is.na(cleanDat))) > LThalfSamples)] #handles case with one bad protein!
removedRownames## character(0)# List of Proteins to be removed:
#635 proteins are removed here!
# remove rows with >=50% missing values (only if there are some rows to be removed)
if (length(na.omit(match(removedRownames, rownames(cleanDat)))) == length(removedRownames) & length(removedRownames) > 0 ) {
cleanDat <- cleanDat[-match(removedRownames, rownames(cleanDat)), ]
} else {
cat("")
}
dim(cleanDat)## [1] 1544 15#Check for genes and samples with too many missing values
gag = goodSamplesGenes(cleanDat, verbose = 3)## Flagging genes and samples with too many missing values...
## ..step 1gag$allOK## [1] TRUEsampleTree = hclust(dist(t(cleanDat)), method = "average")
# Plot sample tree
sizeGrWindow(12, 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 = 15, col = "red")
clust = cutreeStatic(sampleTree, cutHeight = 15, minSize = 10)
table(clust)## clust
## 0
## 15library("doParallel")## Loading required package: foreach## Loading required package: iterators## Loading required package: parallel#stopCluster(clusterLocal) #in case already started.
parallelThreads=12
clusterLocal <- makeCluster(c(rep("localhost",parallelThreads)),type="SOCK")
registerDoParallel(clusterLocal)
allowWGCNAThreads() ## Allowing multi-threading with up to 8 threads.#you want an R squared close to 0.8 or more and a mean k of around a 100
powers <- seq(4,25,by=1) #initial power check -- try to get SFT.R.sq to go > 0.80
sft <- pickSoftThreshold(t(cleanDat),blockSize=nrow(cleanDat)+1000, #always calculate power within a single block (blockSize > # of rows in cleanDat)
powerVector=powers,
corFnc="bicor",networkType="signed")## Warning in (function (x, y = NULL, robustX = TRUE, robustY = TRUE, use =
## "all.obs", : bicor: zero MAD in variable 'x'. Pearson correlation was used for
## individual columns with zero (or missing) MAD.## Warning in (function (x, y = NULL, robustX = TRUE, robustY = TRUE, use =
## "all.obs", : bicor: zero MAD in variable 'x'. Pearson correlation was used for
## individual columns with zero (or missing) MAD.## Warning in (function (x, y = NULL, robustX = TRUE, robustY = TRUE, use =
## "all.obs", : bicor: zero MAD in variable 'x'. Pearson correlation was used for
## individual columns with zero (or missing) MAD.## Warning in (function (x, y = NULL, robustX = TRUE, robustY = TRUE, use =
## "all.obs", : bicor: zero MAD in variable 'x'. Pearson correlation was used for
## individual columns with zero (or missing) MAD.## Warning in (function (x, y = NULL, robustX = TRUE, robustY = TRUE, use =
## "all.obs", : bicor: zero MAD in variable 'x'. Pearson correlation was used for
## individual columns with zero (or missing) MAD.## Warning in (function (x, y = NULL, robustX = TRUE, robustY = TRUE, use =
## "all.obs", : bicor: zero MAD in variable 'x'. Pearson correlation was used for
## individual columns with zero (or missing) MAD.## Warning in (function (x, y = NULL, robustX = TRUE, robustY = TRUE, use =
## "all.obs", : bicor: zero MAD in variable 'x'. Pearson correlation was used for
## individual columns with zero (or missing) MAD.## Warning in (function (x, y = NULL, robustX = TRUE, robustY = TRUE, use =
## "all.obs", : bicor: zero MAD in variable 'x'. Pearson correlation was used for
## individual columns with zero (or missing) MAD.## Warning in (function (x, y = NULL, robustX = TRUE, robustY = TRUE, use =
## "all.obs", : bicor: zero MAD in variable 'y'. Pearson correlation was used for
## individual columns with zero (or missing) MAD.## Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k.
## 1 4 0.598 -3.00 0.833 156.00 150.000 227.0
## 2 5 0.821 -2.78 0.943 102.00 95.300 172.0
## 3 6 0.902 -2.53 0.983 69.10 62.600 135.0
## 4 7 0.932 -2.31 0.985 48.50 42.400 109.0
## 5 8 0.946 -2.15 0.987 35.00 29.400 90.3
## 6 9 0.956 -2.04 0.994 25.90 20.700 76.0
## 7 10 0.957 -1.95 0.992 19.60 14.900 65.0
## 8 11 0.959 -1.88 0.990 15.10 10.900 56.3
## 9 12 0.967 -1.83 0.997 11.90 8.130 49.2
## 10 13 0.969 -1.80 0.994 9.45 6.150 43.5
## 11 14 0.954 -1.78 0.981 7.63 4.690 38.7
## 12 15 0.934 -1.77 0.964 6.23 3.620 34.6
## 13 16 0.942 -1.73 0.969 5.14 2.820 31.2
## 14 17 0.937 -1.72 0.966 4.29 2.220 28.3
## 15 18 0.943 -1.70 0.971 3.61 1.750 25.7
## 16 19 0.931 -1.69 0.960 3.06 1.400 23.5
## 17 20 0.939 -1.67 0.965 2.62 1.130 21.6
## 18 21 0.943 -1.64 0.970 2.25 0.917 19.9
## 19 22 0.947 -1.62 0.972 1.95 0.744 18.4
## 20 23 0.951 -1.60 0.974 1.70 0.613 17.0
## 21 24 0.954 -1.59 0.975 1.49 0.506 15.8
## 22 25 0.959 -1.57 0.980 1.31 0.420 14.7saveRDS(sft, file = 'sft.rds')
###plot
sizeGrWindow(9,5)
par(mfrow= c(1,2))
cex1=0.9
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")
abline(h=0.90, 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")library(ggplot2)
library(RColorBrewer)
#plot initial SFT.R.sq vs. power curve
tableSFT<-sft[[2]]
plot(tableSFT[,1],tableSFT[,2],xlab="Power (Beta)",ylab="SFT R")
softPower = sft$powerEstimate
adjacency = adjacency(t(cleanDat), power = softPower)
TOM = TOMsimilarity(adjacency)## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.dissTOM = 1 - TOM
geneTree = hclust(as.dist(dissTOM), method = "average")
dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM, deepSplit = 2, pamRespectsDendro = FALSE)## ..cutHeight not given, setting it to 0.993 ===> 99% of the (truncated) height range in dendro.
## ..done.moduleColors = labels2colors(dynamicMods)
MEs = moduleEigengenes(t(cleanDat), colors = moduleColors)$eigengenes
MEs = orderMEs(MEs)
nSamples = nrow(t(cleanDat))
# Convert traits to numeric
numericTraits = as.data.frame(lapply(datTraits, function(x) as.numeric(as.factor(x))))
moduleTraitCor = cor(MEs, numericTraits, use = "p")
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)
#Define a function to create plots with r and p values:
plotModuleTrait = function(MEname, trait) {
dat = data.frame(ME = MEs[, MEname], Trait = as.factor(datTraits[, trait]))
cor_value = cor(dat$ME, as.numeric(dat$Trait))
p_value = corPvalueStudent(cor_value, nSamples)
plot_title = paste(MEname, "\nr =", round(cor_value, 2), "p =", signif(p_value, 2))
p = ggplot(dat, aes(x = Trait, y = ME, fill = Trait)) +
geom_boxplot() +
theme_minimal() +
labs(title = plot_title, x = trait, y = "Module Eigengene") +
stat_summary(fun.data = mean_cl_normal, geom = "errorbar", width = 0.2) +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
print(p)
}
# Open a PDF device
pdf("Module_Trait_Plots.pdf", width = 8, height = 6)
# Generate plots for all relevant modules and traits
moduleNames = names(MEs)
traitNames = names(datTraits)
for (module in moduleNames) {
for (trait in traitNames) {
plotModuleTrait(module, trait)
}
}
# Close the PDF device
dev.off()## quartz_off_screen
## 2power<-powerSFT<- 12
getwd()## [1] "/Users/usri/Desktop/Camk2a.WGCNA.8june/JUNE30"rootdir <- "/Users/usri/Desktop/Camk2a.WGCNA.8june/JUNE30/"
minModSize= 13 #for protein; in RNA, this is still 25 #specifies that min number of prots that need to be in a module!
enforceMMS=TRUE
#Read up on blockwise Modules
##making of net modules
net <- blockwiseModules(t(cleanDat),power=power,deepSplit=2,minModuleSize=minModSize,
mergeCutHeight=0.07, TOMDenom="mean", #detectCutHeight=0.9999,
corType="bicor",networkType="signed",pamStage=TRUE,pamRespectsDendro=TRUE,reassignThresh=0.05,
verbose=3,saveTOMs=FALSE,maxBlockSize=40000,nThreads=15)## Calculating module eigengenes block-wise from all genes
## Flagging genes and samples with too many missing values...
## ..step 1
## ..Working on block 1 .
## TOM calculation: adjacency..
## ..will use 15 parallel threads.
## Fraction of slow calculations: 0.000000
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.
## ....clustering..
## ....detecting modules..
## ....calculating module eigengenes..
## ....checking kME in modules..
## ..removing 12 genes from module 1 because their KME is too low.
## ..removing 8 genes from module 2 because their KME is too low.
## ..removing 14 genes from module 3 because their KME is too low.
## ..removing 13 genes from module 4 because their KME is too low.
## ..removing 12 genes from module 5 because their KME is too low.
## ..removing 12 genes from module 6 because their KME is too low.
## ..removing 6 genes from module 7 because their KME is too low.
## ..removing 4 genes from module 8 because their KME is too low.
## ..removing 6 genes from module 9 because their KME is too low.
## ..removing 3 genes from module 10 because their KME is too low.
## ..removing 2 genes from module 12 because their KME is too low.
## ..reassigning 7 genes from module 1 to modules with higher KME.
## ..reassigning 3 genes from module 2 to modules with higher KME.
## ..reassigning 3 genes from module 3 to modules with higher KME.
## ..reassigning 1 genes from module 4 to modules with higher KME.
## ..reassigning 1 genes from module 8 to modules with higher KME.
## ..reassigning 1 genes from module 10 to modules with higher KME.
## ..merging modules that are too close..
## mergeCloseModules: Merging modules whose distance is less than 0.07
## Calculating new MEs...saveRDS(net, file = 'net.rds')
nModules<-length(table(net$colors))-1
modules<-cbind(colnames(as.matrix(table(net$colors))),table(net$colors))
orderedModules<-cbind(Mnum=paste("M",seq(1:nModules),sep=""),Color=labels2colors(c(1:nModules)))
modules<-modules[match(as.character(orderedModules[,2]),rownames(modules)),]
MEs<-tmpMEs<-data.frame()
MEList = moduleEigengenes(t(cleanDat), colors = net$colors)
MEs = orderMEs(MEList$eigengenes)
net$MEs <- MEs
colnames(MEs)<-gsub("ME","",colnames(MEs)) #let's be consistent in case prefix was added, remove it.
rownames(MEs)<-rownames(numericMeta)
tmpMEs <- MEs #net$MEs
colnames(tmpMEs) <- paste("ME",colnames(MEs),sep="")
MEs[,"grey"] <- NULL
tmpMEs[,"MEgrey"] <- NULL
kMEdat <- signedKME(t(cleanDat), tmpMEs, corFnc="bicor")## Warning in bicor(datExpr, datME, , use = "p"): bicor: zero MAD in variable 'x'.
## Pearson correlation was used for individual columns with zero (or missing) MAD.#saveRDS(kMEdat, file = 'kMEdat.rds')
#The protins that do not fall into the kme calc fall into module gray and are excluded!FileBaseName <- 'Camk2a.30june'
rootdir <- "/Users/usri/Desktop/Camk2a.WGCNA.8june/JUNE30"
outputtabs<-outputfigs<-rootdir
#(re)Calculate MEs (no grey)
MEs<-tmpMEs<-data.frame()
MEList = moduleEigengenes(t(cleanDat), colors = net$colors)
MEs = orderMEs(MEList$eigengenes)
colnames(MEs)<-gsub("ME","",colnames(MEs))
tmpMEs <- MEs #net$MEs
colnames(tmpMEs) <- paste("ME",colnames(MEs),sep="")
MEs[,"grey"] <- NULL
tmpMEs[,"MEgrey"] <- NULL
#specify filename
baseName=FileBaseName
MEs <- net$MEs
#(re)Calculate MEs (no grey)
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.
tmpMEs <- MEs #net$MEs
colnames(tmpMEs) <- paste("ME",colnames(MEs),sep="")
MEs[,"grey"] <- NULL
tmpMEs[,"MEgrey"] <- NULL
#specify filename
baseName=FileBaseName #inputExprMat #PDF filename base
#check these settings
#######################################################
PPIedges=FALSE #TRUE will be a lot slower...
vertexsize=16 #8 for regular, 16 for large balls
species="human" #current option "mouse" will convert bioGRID to mouse symbols before drawing PPI edges.
#CAIRO=TRUE #Open/Write PDF using CairoPDF or pdf functions; now both file versions are created
########################################################
#(re)Establish M# for module colors table
nModules<-length(table(net$colors))-1
modules<-cbind(colnames(as.matrix(table(net$colors))),table(net$colors))
orderedModules<-cbind(Mnum=paste("M",seq(1:nModules),sep=""),Color=labels2colors(c(1:nModules)))
#load interactome [full bioGRID 3.4 (01/25/2016)]
#bioGrid <- read.table(file=paste0(datadir,"PPIs/BIOGRID-Homo_sapiens-3.4.133.SIMPLEsymbols_editpad.txt"),sep="\t",header=TRUE)
#------
if (species=="mouse") {
library(biomaRt)
human = useMart("ENSEMBL_MART_ENSEMBL",dataset="hsapiens_gene_ensembl",host="https://www.ensembl.org")
mouse = useMart("ENSEMBL_MART_ENSEMBL",dataset="mmusculus_gene_ensembl",host="https://www.ensembl.org")
## If converting bioGrid from one species to another
# genelist.convert<-getLDS(attributes=c("hgnc_symbol"), filters="hgnc_symbol", values=unique(c(as.vector(bioGrid[,"From"]),as.vector(bioGrid[,"To"]))), mart=human, attributesL=c("mgi_symbol","external_gene_name"),martL = mouse)
#
# bioGrid.mouse.convert<-cbind(genelist.convert[match(bioGrid[,1],genelist.convert[,1]),3],genelist.convert[match(bioGrid[,2],genelist.convert[,1]),3])
# bioGrid.human<-bioGrid
# bioGrid<-bioGrid.mouse.convert
# dim(bioGrid)
# bioGrid<-as.data.frame(bioGrid[-which(is.na(bioGrid)),])
# colnames(bioGrid)<-colnames(bioGrid.human)
# dim(bioGrid)
}
if (PPIedges==FALSE) { bioGrid<-data.frame(From=c(1),To=c(1)) } #if you want to skip bolding edges for PPIs from BioGrid
#Symbol to add to iGraphs that don't include it, if checking for protein-protein interaction edges (add 4, one for each corner; can edit output later to remove unwanted ones.)
symbols2Add=c() #c("PiB","APP","SGIP1","APOE")
RNAbindingGOlist<- c()
#Get KMEs
# KME.vis=signedKME(cleanDat, kMEdat,corFnc="bicor",outputColumnName = "");
KME.vis=kMEdat[,] # -ncol(kMEdat) #minus grey column if calculating signedKME on the fly
annot=as.data.frame(do.call("rbind",strsplit(as.character(rownames(cleanDat)),"[|]")))
KME.vis$Symbol=annot$V1
geneInfo.vis=as.data.frame(cbind(KME.vis$Symbol,net$colors, KME.vis))
geneInfo.vis=geneInfo.vis[,-ncol(geneInfo.vis)] # check if last column is Ensembl gene id
colnames(geneInfo.vis)[1]= "Symbol"
colnames(geneInfo.vis)[2]= "Module.Color"
colnames(geneInfo.vis)<-gsub("kMEME","kME",colnames(geneInfo.vis))
library(igraph);## Warning: package 'igraph' was built under R version 4.3.2##
## Attaching package: 'igraph'## The following object is masked from 'package:IRanges':
##
## union## The following object is masked from 'package:S4Vectors':
##
## union## The following objects are masked from 'package:BiocGenerics':
##
## normalize, path, union## The following objects are masked from 'package:dplyr':
##
## as_data_frame, groups, union## The following object is masked from 'package:tidyr':
##
## crossing## The following objects are masked from 'package:stats':
##
## decompose, spectrum## The following object is masked from 'package:base':
##
## union#iGRAPHs Connectivity Hub circular Plot
library(RColorBrewer);
library(WGCNA);
library(Cairo);
softPower = power
adjacency = adjacency(cleanDat, power=softPower, type="signed",corFnc="bicor")
TOM = TOMsimilarity(adjacency)## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.TOM.matrix = as.matrix(TOM);
#Get the top connected genes in the module
uniquemodcolors = gsub("kME","",gsub("kMEME","",colnames(kMEdat[,])));
for (CAIRO in c(TRUE,FALSE)) {
if (CAIRO) { CairoPDF(file=paste0(outputfigs,"/6.iGraph.camk2_Modules-",FileBaseName,"-CAIRO-LargeNodes.pdf"),width=16,height=12) } else {
pdf(paste0(outputfigs,"/6.iGraph.camk2a_Modules-",FileBaseName,"-nonCAIRO-LargeNodes.pdf"),height=9,width=10)
}
# for (i in 1:length(sigmodcolors)) {
# mod=sigmodcolors[i];
# numgenesingraph = 50;
# numconnections2keep = 1500;
for (mod in uniquemodcolors) {
#mod="darkslateblue"
numgenesingraph = 100;
numconnections2keep = 700;
cat('module:',mod,'\n');
geneInfo.vis=geneInfo.vis[geneInfo.vis$Symbol!="NA",]
geneInfo.vis=geneInfo.vis[geneInfo.vis$Symbol!="",]
colind = which(colnames(geneInfo.vis)== paste("kME",mod,sep=""));
rowind = which(geneInfo.vis[,2]==mod);
cat(' ',length(rowind),'probes in module\n');
submatrix = geneInfo.vis[rowind,];
orderind = order(submatrix[,colind],decreasing=TRUE);
if (length(rowind) < numgenesingraph) {
numgenesingraph = length(rowind);
numconnections2keep = numgenesingraph/2 * (numgenesingraph/6 - 1); #added /2 and /6 9/14/2015
}
if (length(rowind)<10) { innercircleNum=length(rowind) } else { innercircleNum=10 }
cat('Making network graphs, using top',numgenesingraph,'probes and',numconnections2keep,'connections of TOM\n');
submatrix = submatrix[orderind[1:numgenesingraph],];
#Identify the columns in the TOM that correspond to these hub probes
matchind = match(submatrix$Symbol,annot$Symbol);
reducedTOM = TOM.matrix[matchind,matchind];
orderind = order(reducedTOM,decreasing=TRUE);
connections2keep = orderind[1:numconnections2keep];
reducedTOM = matrix(0,nrow(reducedTOM),ncol(reducedTOM));
reducedTOM[connections2keep] = 1;
g0 <- graph.adjacency(as.matrix(reducedTOM[1:innercircleNum,1:innercircleNum]),mode="undirected",weighted=TRUE,diag=FALSE)
layoutMata <- layout.circle(g0)
if (ncol(reducedTOM) < 51 & ncol(reducedTOM) > 10) {
g0 <- graph.adjacency(as.matrix(reducedTOM[11:ncol(reducedTOM),11:ncol(reducedTOM)]),mode="undirected",weighted=TRUE,diag=FALSE)
layoutMatb <- layout.circle(g0)
g1 <- graph.adjacency(as.matrix(reducedTOM),mode="undirected",weighted=TRUE,diag=FALSE)
layoutMat <- rbind(layoutMata*0.25,layoutMatb*0.75)
} else { if (ncol(reducedTOM) > 10) { #****
g0 <- graph.adjacency(as.matrix(reducedTOM[11:50,11:50]),mode="undirected",weighted=TRUE,diag=FALSE)
layoutMatb <- layout.circle(g0)
g0 <- graph.adjacency(as.matrix(reducedTOM[51:ncol(reducedTOM),51:ncol(reducedTOM)]),mode="undirected",weighted=TRUE,diag=FALSE)
layoutMatc <- layout.circle(g0)
g1 <- graph.adjacency(as.matrix(reducedTOM),mode="undirected",weighted=TRUE,diag=FALSE)
layoutMat <- rbind(layoutMata*0.25,layoutMatb*0.75, layoutMatc)
}
else { layoutMat <- layoutMata }
}
#Add symbols2Add
iter=as.numeric(0)
if (length(symbols2Add)==4) {
for (symbol2Add in symbols2Add) {
iter=iter+1
if (iter==1) { position=cbind(-0.7071068,-0.7071068) }
if (iter==2) { position=rbind(position,cbind(0.7071068,-0.7071068)) }
if (iter==3) { position=rbind(position,cbind(-0.7071068,0.7071068)) }
if (iter==4) { position=rbind(position,cbind(0.7071068,0.7071068)) }
#ADD APP to graph g3 (to be used if APP not already in graph)
g2 <- vertex(1) #,size=40,color="green", label=symbol2Add
#layoutMat <- rbind(layoutMata*0.25,rbind(layoutMatb*0.75,cbind(-0.7071068,-0.7071068)))
g3 <- g1+g2
V(g3)$name <- c(1:(numgenesingraph+1))
#WORKS# plot(g3,edge.color="grey",vertex.color=c(rep(mod,nrow(layoutMat)-1),"green"),vertex.label=as.character(c(submatrix$Symbol,symbol2Add)),vertex.label.cex=1.1,vertex.label.dist=0.45,vertex.label.degree=-pi/4,vertex.label.color="black",layout= layoutMat,vertex.size=c(submatrix[,colind]^2*16,15),main=paste(mod,"module"))
g1<-g3
numgenesingraph=numgenesingraph+1
}
}
#x***moved out of loop
if (length(symbols2Add)==4) { if (ncol(reducedTOM) < 51) { layoutMat <- rbind(layoutMata*0.25,rbind(layoutMatb*0.75,position)) } else { layoutMat <- rbind(layoutMata*0.25,layoutMatb*0.75, rbind(layoutMatc,position*1.33)) }
} else { if (ncol(reducedTOM) < 51 & ncol(reducedTOM) > 10) { layoutMat <- rbind(layoutMata*0.25,layoutMatb*0.75) } else { if(ncol(reducedTOM) > 10) { layoutMat <- rbind(layoutMata*0.25,layoutMatb*0.75, layoutMatc) } } #do not add 4 node positions if symbols2Add is empty (length 0)
}
symbolList<-c(as.character(submatrix$Symbol),symbols2Add)
if (length(symbols2Add)==4) { vertexColorVec<-c(rep(mod,numgenesingraph)-4,"green","darkgreen","steelblue","darkslateblue"); vertexSizeMat<-c(submatrix[,colind]^2*vertexsize,rep(15,4)); } else { vertexColorVec<-rep(mod,numgenesingraph); vertexSizeMat<-submatrix[,colind]^2*vertexsize; }
#FIND EDGES OVERLAPPING WITH BIOGRID PAIRS
listboldEdges<-matrix(ncol=2,nrow=0)
if(nrow(bioGrid)>0) {
for(i in 1:numgenesingraph) {
for(j in i:numgenesingraph) {
if(!(length(which(bioGrid$From==symbolList[i]&bioGrid$To==symbolList[j]))+length(which(bioGrid$From==symbolList[j]&bioGrid$To==symbolList[i])))==0) { listboldEdges <- rbind(listboldEdges,c(i,j)) }
} }
##Remove self-loops
if (PPIedges) { if (length(nrow(listboldEdges))>0 & nrow(listboldEdges)>0) {
for(i in nrow(listboldEdges):1) {
if (i<=nrow(listboldEdges)) {
if(listboldEdges[i,1]==listboldEdges[i,2]) {
listboldEdges <-listboldEdges[-i,]
if (!length(nrow(listboldEdges))>0) {
cat('NO PHYSICAL INTERACTIONS FOUND FOR THIS MODULE.\n')
break #HANDLE NO BOLD EDGES CASE
}
}
}
}
}
}
}
if (is.vector(listboldEdges)) { listboldEdges<-t(data.frame(listboldEdges)) }
newgraph <- g1 %>%
# delete_edges(listboldEdges) %>%
set_edge_attr("color",value="lightgrey") %>%
set_edge_attr("width",value=1) %>%
set_edge_attr("curved",value=0)
if (length(symbols2Add)==4) {
if (is.vector(listboldEdges)) {
cat('ONLY 1 SINGLE INTERACTION FOUND.\n')
newgraph <- newgraph %>%
add_edges(c(listboldEdges[1],listboldEdges[2]),color="steelblue",width=2,curved=0)
} else {
if (dim(listboldEdges)[1]>0) {
for(k in 1:nrow(listboldEdges)) {
if((length(which(submatrix$Symbol==symbols2Add[1]))==0)&(listboldEdges[k,1]==(numgenesingraph-3)|listboldEdges[k,2]==(numgenesingraph-3))) {
newgraph <- newgraph %>%
add_edges(c(listboldEdges[k,1],listboldEdges[k,2]),color="#33BB33",width=3,curved=0)
}
}
for(k in 1:nrow(listboldEdges)) {
if((length(which(submatrix$Symbol==symbols2Add[2]))==0)&(listboldEdges[k,1]==(numgenesingraph-2)|listboldEdges[k,2]==(numgenesingraph-2))) {
newgraph <- newgraph %>%
add_edges(c(listboldEdges[k,1],listboldEdges[k,2]),color="#338833",width=3,curved=0)
}
}
for(k in 1:nrow(listboldEdges)) {
if((length(which(submatrix$Symbol==symbols2Add[3]))==0)&(listboldEdges[k,1]==(numgenesingraph-1)|listboldEdges[k,2]==(numgenesingraph-1))) {
newgraph <- newgraph %>%
add_edges(c(listboldEdges[k,1],listboldEdges[k,2]),color="steelblue",width=3,curved=0)
}
}
for(k in 1:nrow(listboldEdges)) {
if((length(which(submatrix$Symbol==symbols2Add[4]))==0)&(listboldEdges[k,1]==numgenesingraph|listboldEdges[k,2]==numgenesingraph)) {
newgraph <- newgraph %>%
add_edges(c(listboldEdges[k,1],listboldEdges[k,2]),color="darkslateblue",width=3,curved=0)
}
}
}
}
} else {
if (dim(listboldEdges)[1]>0) {
for(k in 1:nrow(listboldEdges)) {
newgraph <- newgraph %>%
add_edges(c(listboldEdges[k,1],listboldEdges[k,2]),color="#33BB33",width=3,curved=0)
}
}
}
highlightNodes<-match(RNAbindingGOlist,symbolList) #ADmagmaList$Gene.Symbol,symbolList
highlightNodes<-sort(na.omit(highlightNodes))
if(length(highlightNodes)>0) { if(mod=="yellow") {vertexColorVec[highlightNodes]<-"cyan" } else {vertexColorVec[highlightNodes]<-"yellow" } }
plot(newgraph,vertex.color=vertexColorVec,vertex.label=as.character(symbolList),vertex.label.cex=1.1,vertex.label.dist=0.45,vertex.label.degree=-pi/4,vertex.label.color="black",layout=layoutMat,vertex.size=vertexSizeMat,main=paste0(orderedModules[which(orderedModules[,"Color"]==mod),"Mnum"]," ",mod," module"))
}
dev.off();
}## module: tan
## 21 probes in module
## Making network graphs, using top 21 probes and 26.25 connections of TOM## Warning: `graph.adjacency()` was deprecated in igraph 2.0.0.
## ℹ Please use `graph_from_adjacency_matrix()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## module: turquoise
## 210 probes in module
## Making network graphs, using top 100 probes and 700 connections of TOM## module: yellow
## 126 probes in module
## Making network graphs, using top 100 probes and 700 connections of TOM## module: black
## 58 probes in module
## Making network graphs, using top 58 probes and 251.3333 connections of TOM## module: purple
## 49 probes in module
## Making network graphs, using top 49 probes and 175.5833 connections of TOM## module: greenyellow
## 38 probes in module
## Making network graphs, using top 38 probes and 101.3333 connections of TOM## module: pink
## 57 probes in module
## Making network graphs, using top 57 probes and 242.25 connections of TOM## module: green
## 109 probes in module
## Making network graphs, using top 100 probes and 700 connections of TOM## module: blue
## 206 probes in module
## Making network graphs, using top 100 probes and 700 connections of TOM## module: red
## 108 probes in module
## Making network graphs, using top 100 probes and 700 connections of TOM## module: salmon
## 20 probes in module
## Making network graphs, using top 20 probes and 23.33333 connections of TOM## module: brown
## 130 probes in module
## Making network graphs, using top 100 probes and 700 connections of TOM## module: magenta
## 54 probes in module
## Making network graphs, using top 54 probes and 216 connections of TOM## module: tan
## 21 probes in module
## Making network graphs, using top 21 probes and 26.25 connections of TOM## module: turquoise
## 210 probes in module
## Making network graphs, using top 100 probes and 700 connections of TOM## module: yellow
## 126 probes in module
## Making network graphs, using top 100 probes and 700 connections of TOM## module: black
## 58 probes in module
## Making network graphs, using top 58 probes and 251.3333 connections of TOM## module: purple
## 49 probes in module
## Making network graphs, using top 49 probes and 175.5833 connections of TOM## module: greenyellow
## 38 probes in module
## Making network graphs, using top 38 probes and 101.3333 connections of TOM## module: pink
## 57 probes in module
## Making network graphs, using top 57 probes and 242.25 connections of TOM## module: green
## 109 probes in module
## Making network graphs, using top 100 probes and 700 connections of TOM## module: blue
## 206 probes in module
## Making network graphs, using top 100 probes and 700 connections of TOM## module: red
## 108 probes in module
## Making network graphs, using top 100 probes and 700 connections of TOM## module: salmon
## 20 probes in module
## Making network graphs, using top 20 probes and 23.33333 connections of TOM## module: brown
## 130 probes in module
## Making network graphs, using top 100 probes and 700 connections of TOM## module: magenta
## 54 probes in module
## Making network graphs, using top 54 probes and 216 connections of TOM#net <- readRDS("~/Desktop/Camk2a.WGCNA.8june/JUNE30/net.rds")
#MEs <- net$MEs
#(re-Calculate MEs(no grey)
#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.
#tmpMEs <- MEs #net$MEs
#colnames(tmpMEs) <- paste("ME",colnames(MEs),sep="")
#MEs[,"grey"] <- NULL
#tmpMEs[,"MEgrey"] <- NULL
MEs.no.gray <- readRDS("~/Desktop/Camk2a.WGCNA.8june/JUNE30/MEs.no.gray.rds")
numericMeta <- readRDS("~/Desktop/Camk2a.WGCNA.8june/JUNE30/numericMeta.rds")
MEs <- MEs.no.gray
regvars <- data.frame(as.factor(numericMeta[,"Genotype"]))
colnames(regvars) <- c("Genotype")
lm <- lm(data.matrix(MEs)~Genotype,data=regvars)
pvec <- rep(NA,ncol(MEs))
for (i in 1:ncol(MEs)) {
f <- summary(lm)[[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)
Groupingx <- numericMeta$Genotype
metdat <- data.frame(Status=Groupingx)
#toplot <- MEs
#colnames(toplot) <- colnames(MEs)
#rownames(toplot) <- rownames(MEs)
#toplot <- t(toplot) #may have to transpose sometimes: t(toplot)
#rownames(toplot) <- paste(orderedModules[match(colnames(MEs),orderedModules[,2]),1]," ",rownames(toplot),sep="")
#I save this, so here I directly load from source
toplot <- readRDS("~/Desktop/Camk2a.WGCNA.8june/JUNE30/toplot.rds")
head(toplot, 4)## X.106P2fractionIP X.120P2fractionIP X.117P2fractionIP
## tan 0.002286292 -0.12085997 -0.1836851
## turquoise -0.140061765 -0.03180121 -0.1086136
## yellow -0.279556163 -0.32547276 -0.4412361
## greenyellow -0.002570635 0.40846040 0.4238847
## X.123P2fractionIP X.137P2fractionIP X.131P2fractionIP
## tan -0.201417732 -0.41213175 0.3594462
## turquoise -0.009744013 -0.38341634 -0.3987381
## yellow -0.342866509 -0.08059125 -0.1038111
## greenyellow 0.412058843 -0.26524084 -0.3604647
## X.143P2fractionIP X.133P2fractionIP X.106HomogenateIP
## tan -0.02059604 -0.3961085 0.23177639
## turquoise -0.30103479 -0.3684375 0.23685035
## yellow -0.06961930 -0.1323982 0.23919658
## greenyellow -0.15330009 -0.3176283 0.04371573
## X.120HomogenateIP X.117HomogenateIP X.123HomogenateIP
## tan -0.07633985 -0.1957735 0.04341582
## turquoise 0.23516251 0.2391938 0.23858100
## yellow 0.25794282 0.1777519 0.27887421
## greenyellow -0.13874270 -0.1538161 0.09215113
## X.131HomogenateIP X.143HomogenateIP X.133HomogenateIP
## tan 0.3166204 0.4702886 0.1830788
## turquoise 0.2743220 0.2710859 0.2466517
## yellow 0.2348201 0.3353184 0.2516472
## greenyellow -0.1167498 0.2688664 -0.1406241par(mfrow=c(3,3))
par(mar=c(4.5,6,4.5,1.5))
#pdf(file=paste0(outputfigs,"/BoxPlots_",FileBaseName,".pdf"),width=18,height=11.25)
# Replace this with your actual data or method of creating orderedModules
orderedModules <- data.frame(
ModuleName = colnames(MEs),
Color = colnames(MEs))
colnames(orderedModules) <- c("ModuleName", "Color")
orderedModules## ModuleName Color
## 1 tan tan
## 2 turquoise turquoise
## 3 yellow yellow
## 4 greenyellow greenyellow
## 5 pink pink
## 6 magenta magenta
## 7 purple purple
## 8 green green
## 9 blue blue
## 10 red red
## 11 brown brown
## 12 black black
## 13 salmon salmonlibrary(beeswarm)
library(gplots)## Warning: package 'gplots' was built under R version 4.3.2##
## Attaching package: 'gplots'## The following object is masked from 'package:IRanges':
##
## space## The following object is masked from 'package:S4Vectors':
##
## space## The following object is masked from 'package:stats':
##
## lowessfor (i in 1:(nrow(toplot))) { # grey already excluded, no -1
titlecolor<-if(signif(pvec,2)[i] <0.05) { "red" } else { "black" }
boxplot(toplot[i,]~factor(numericMeta$Genotype),col=colnames(MEs)[i],ylab="Eigenprotein Value",main=paste0(orderedModules[match(colnames(MEs)[i],orderedModules[,2]),1]," ",colnames(MEs)[i],"\nANOVA p = ",signif(pvec,2)[i]),xlab="Genotype",las=1,col.main=titlecolor) #no outliers: ,outline=FALSE)
transcol=paste0(col2hex(colnames(MEs)[i]),"99")
beeswarm(toplot[i,]~factor(numericMeta$Genotype),method="swarm",add=TRUE,corralWidth=0.5,vertical=TRUE,pch=21,bg=transcol,col="black",cex=0.8,corral="gutter") }
