RNAseq|WGCNA-组学数据黏合剂，代码实战-一（尽）文（力）解决文献中常见的可视化图

library(tidyverse)
library(openxlsx)
library(data.table)
#BiocManager::install("impute")
library(WGCNA)

load("RNAseq.SKCM_WGCNA.RData")
head(cibersort_in,2)
expr2[1:4,1:4]

## 选取top 5000 mad genes
WGCNA_matrix = t(expr2[order(apply(expr2,1,mad), decreasing = T)[1:5000],]) 
datExpr <- WGCNA_matrix  

### beta value#########
powers = c(c(1:10), seq(from = 12, to=30, by=2))
sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)

par(mfrow = c(1,2))
cex1 = 0.9; #可自行调整
# Scale-free topology fit index as a function of the soft-thresholding power
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.90,col="red")
# Mean connectivity as a function of the soft-thresholding power
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")

sft$powerEstimate

net = blockwiseModules(
  datExpr,
  power = sft$powerEstimate,
  #power = sft$powerEstimate,
  maxBlockSize = 6000,
  TOMType = "unsigned", minModuleSize = 30,
  reassignThreshold = 0, mergeCutHeight = 0.25,
  numericLabels = TRUE, pamRespectsDendro = FALSE,
  saveTOMs = TRUE,
  saveTOMFileBase = "AS-green-FPKM-TOM",
  verbose = 3
)
table(net$colors) 
####模块可视化####
mergedColors = labels2colors(net$colors)
table(mergedColors)
# 可视化
plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
                    "Module colors",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)

#随机选取部分基因作图
nGenes = ncol(datExpr)
nSamples = nrow(datExpr)
nSelect = 500
# 设置随机种子，方便复现结果
set.seed(1234);
select = sample(nGenes, size = nSelect);
selectTOM = dissTOM[select, select];
# There’s no simple way of restricting a clustering tree to a subset of genes, so we must re-cluster.
selectTree = hclust(as.dist(selectTOM), method = "average")
selectColors = moduleColors[select];
# Open a graphical window
sizeGrWindow(9,9)

plotDiss = selectTOM^7;
diag(plotDiss) = NA;

TOMplot(plotDiss, selectTree, selectColors,terrainColors = T, 
        main = "Network heatmap plot, selected genes")

library(gplots)
TOMplot(plotDiss, selectTree, selectColors, main = "Network heatmap plot, selected genes", 
        col=gplots::colorpanel(250,'red',"orange",'lemonchiffon'))

datTraits <- cibersort_in %>% 
  select(sample:Neutrophils ) %>%
  column_to_rownames("sample")
nSamples = nrow(datExpr)
sampleNames = rownames(datExpr)
traitRows = match(datTraits$sample, rownames(datTraits) )  

moduleColors <- labels2colors(net$colors)
# Recalculate MEs with color labels
MEs0 = moduleEigengenes(datExpr, moduleColors)$eigengenes
MEs = orderMEs(MEs0);
##计算相关系数和P值
moduleTraitCor = cor(MEs, datTraits , use = "p");
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)

sizeGrWindow(10,6)
# 图中展示P值
textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
                   signif(moduleTraitPvalue, 1), ")", sep = "");
dim(textMatrix) = dim(moduleTraitCor)
#通过par设置图形边界
par(mar = c(10, 10, 3, 3));
# 热图的形式展示结果
labeledHeatmap(Matrix = moduleTraitCor,
               xLabels = colnames(datTraits),
               yLabels = names(MEs),
               ySymbols = names(MEs),
               colorLabels = FALSE,
               colors = blueWhiteRed(50),
               textMatrix = textMatrix,
               setStdMargins = FALSE,
               cex.text = 0.5,
               zlim = c(-1,1),
               main = paste("Module-trait relationships"))

########## 显著模块柱形图#######
moduleColors = labels2colors(net$colors)
y <- datTraits[, "T cells CD8"]
GS <- as.numeric(cor(y ,datExpr, use="p"))
GeneSignificance <-  abs(GS)
ModuleSignificance <- tapply(
  GeneSignificance,
  moduleColors, mean, na.rm=T)
plotModuleSignificance(GeneSignificance, moduleColors,ylim = c(0,0.4))

# Recalculate module eigengenes
MEs = moduleEigengenes(datExpr, moduleColors)$eigengenes
CD8_T = as.data.frame(datTraits[,4]);
names(CD8_T) = "CD8_T"
modNames = substring(names(MEs), 3)
# Add the weight to existing module eigengenes
MET = orderMEs(cbind(MEs, CD8_T))
# Plot the relationships among the eigengenes and the trait
sizeGrWindow(5,7.5);
par(cex = 0.9)
plotEigengeneNetworks(MET, "", marDendro = c(0,4,1,2), marHeatmap = c(3,4,1,2), cex.lab = 0.8, xLabelsAngle
                      = 90)

pheno = "CD8_T"
# 获取关注的列
module_column = match(module, modNames)
pheno_column = match(pheno,colnames(datTraits))
# 获取模块内的基因
moduleGenes = moduleColors == module

sizeGrWindow(7, 7);
par(mfrow = c(1,1));
verboseScatterplot(abs(geneModuleMembership[moduleGenes, 
                                            module_column]),
                   abs(geneTraitSignificance[moduleGenes, 1]), 
                   xlab = paste("Module Membership in", module, "module"), 
                   ylab = paste("Gene significance for", pheno), 
                   main = paste("Module membership  vs gene significance\n"), 
                   cex.main = 1.2, cex.lab = 1.2, cex.axis = 
                     1.2, col = module)

module = "blue";
# Select module probes
geneID= colnames(datExpr) ## 我们例子里面的probe就是基因名
inModule = (moduleColors==module);
module_gene = geneID[inModule];
write.csv(module_gene,"module_gene.csv")

hub<- abs(geneModuleMembership$MMblue)>0.8 & abs(geneTraitSignificance)>0.2
hub<-as.data.frame(dimnames(data.frame(datExpr))[[2]][hub])
head(hub,2)
#  dimnames(data.frame(datExpr))[[2]][hub]
#1                                     IGJ
#2                                   CXCL9

save(module_gene,hub ,MEs,geneModuleMembership,geneTraitSignificance,
     file = "SKCM.WGCNA_MM_GS.RData")