院士团队的WGCNA挖掘文章修改成为癌症转移与否关键模块

rm(list = ls()) 
load( file =  '../step1-deg/DEG_deseq2.Rdata' ) 
colnames(DEG_deseq2) 
nrDEG=DEG_deseq2[,c("log2FoldChange",   "padj")] 
head(nrDEG) 
colnames(nrDEG)=c('logFC','P.Value')

# 凡是阈值，都是可以自定义 
logFC_t <- 2.5
pvalue_t <- 0.001
gene_up= rownames( nrDEG[with(nrDEG,
                              logFC > logFC_t & P.Value < pvalue_t),])
gene_down=rownames( nrDEG[with(nrDEG,
                               logFC < -logFC_t & P.Value < pvalue_t),]) 
length(gene_up);length(gene_down)
head(gene_up);head(gene_down)

load(file = '../step1-deg/symbol_matrix.Rdata')  
symbol_matrix[1:4,1:4] ## 基因名字的样品，矩阵
dat = log2(edgeR::cpm(symbol_matrix)+1)
# 这个归一化方法不太推荐
library(limma)
library(edgeR)
library(DESeq2) 
dge <- DGEList(counts=symbol_matrix)
dge <- calcNormFactors(dge)
logCPM <- cpm(dge, log=TRUE, prior.count=3)
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(symbol_matrix)
design
v <- voom(dge,design,plot=TRUE, normalize="quantile")
dat = v$E

dat[1:4,1:4]
boxplot(dat,las=2)
table(group_list)

kp = group_list=='case'
# 矩阵，样品过滤，基因过滤
dat = dat[c(gene_up,gene_down),kp]
dim(dat)

#1.确定软阈值----
library(WGCNA)
dat0 <- t(dat)

powers = c(c(1:10), seq(from = 12, to=20, by=2))
# Call the network topology analysis function
sft = pickSoftThreshold(dat0, 
                        powerVector = powers, 
                        verbose = 5)
po <- sft$powerEstimate;po

# Plot the results:
sizeGrWindow(9, 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")
 

#2.一步构建网络----
# 报错：https://blog.csdn.net/liyunfan00/article/details/91686840
cor <- WGCNA::cor
net = blockwiseModules(dat0, power = po,
                       TOMType = "unsigned", minModuleSize = 30,
                       reassignThreshold = 0, mergeCutHeight = 0.25,
                       numericLabels = TRUE, pamRespectsDendro = FALSE,
                       saveTOMs = F, 
                       verbose = 3) #一步构建网络
#改回去
cor<-stats::cor

class(net)
names(net)
table(net$colors)

#(1)保存net相关信息
moduleLabels = net$colors
moduleColors = labels2colors(net$colors)
MEs = net$MEs;
geneTree = net$dendrograms[[1]];
save(net,MEs, moduleLabels, moduleColors, geneTree, 
     file = "Step2networkConstruction-auto.RData")

  0   1   2   3   4   5   6   7   8   9 
509 993 515 341 290 272 209 197 160 122 

#1.Modules-traits relationships----
# Define numbers of genes and samples
dat <- t(dat)
nGenes = ncol(dat)
nSamples = nrow(dat)

# Recalculate MEs with color labels
MEs0 = moduleEigengenes(dat, moduleColors)$eigengenes
MEs0[1:6,1:6]

MEs = orderMEs(MEs0)
MEs[1:6,1:6]

phe = data.table::fread('TCGA-CHOL.GDC_phenotype.tsv.gz',data.table = F)
colnames(phe)
phe=phe[match(rownames(dat),phe$submitter_id.samples),]
# [41] "pathologic_M"                             
# [42] "pathologic_N"                             
# [43] "pathologic_T" 
group_list = phe$pathologic_M
table(group_list)
group_list=ifelse(group_list=='M0','non','M')
group_list=as.factor(group_list)
# 肿瘤，tmn，stage。。。。

design=model.matrix(~0+ group_list)
colnames(design)=levels(group_list)
moduleTraitCor = cor(MEs, design, use = "p")
head(moduleTraitCor)

moduleTraitPvalue = corPvalueStudent(moduleTraitCor, 
                                     nSamples)
head(moduleTraitPvalue)

#可视化
sizeGrWindow(10,6)

# Will display correlations and their p-values
textMatrix =  paste(signif(moduleTraitCor, 2), 
                    "\n(",signif(moduleTraitPvalue, 1),
                    ")", 
                    sep = "");
dim(textMatrix) = dim(moduleTraitCor)
textMatrix[1:6,1:3]

par(mar = c(3, 12, 3, 1));
# Display the correlation values within a heatmap plot
design1 <- as.data.frame(design)
labeledHeatmap(Matrix = moduleTraitCor,
               xLabels = names(design1),
               yLabels = names(MEs),
               ySymbols = names(MEs),
               colorLabels = FALSE,
               colors = blueWhiteRed(50),
               textMatrix = textMatrix,
               setStdMargins = FALSE,
               cex.text = 1,
               zlim = c(-1,1),
               main = paste("Module-trait relationships"))

rm(list = ls())  ## 魔幻操作，一键清空~
options(stringsAsFactors = F)
library(WGCNA)
load( file = "Step2networkConstruction-auto.RData")

table(moduleColors)
group_g=data.frame(gene=names(net$colors),
                   group=moduleColors)
head(group_g)
#save(group_g,file='wgcna_group_g.Rdata') 
library(clusterProfiler)
# Convert gene ID into entrez genes
head(group_g)
tmp <- bitr(group_g$gene, fromType="SYMBOL", 
            toType="ENTREZID", 
            OrgDb="org.Hs.eg.db")
de_gene_clusters=merge(tmp,group_g,by.x='SYMBOL',by.y='gene')
table(de_gene_clusters$group)
head(de_gene_clusters)

gcSample <- split(de_gene_clusters$ENTREZID, 
                               de_gene_clusters$group) 
xx <- compareCluster(gcSample, fun="enrichKEGG",
                     organism="hsa", pvalueCutoff=0.05)
dotplot(xx)