WGCNA

参考

http://www.bio-info-trainee.com/2535.html
WGCNA分析，简单全面的最新教程
https://www.jianshu.com/p/f0409a045dab

image.png
image.png
image.png

library(WGCNA)

options(stringsAsFactors = FALSE)
# 打开多线程
enableWGCNAThreads()
load(file = './Rdata/step2_id_transed.Rdata')

# 常规表达矩阵，log2转换后或
# Deseq2的varianceStabilizingTransformation转换的数据
# 如果有批次效应，需要事先移除，可使用removeBatchEffect
# 如果有系统偏移(可用boxplot查看基因表达分布是否一致)，
# 需要quantile normalization

exprSet<- log2(new_exprSet+1)
pdata<- pdata[1:20,]
exprSet_test<- exprSet[,1:20]
# 官方推荐 "signed" 或 "signed hybrid"
# 为与原文档一致，故未修改 
type = "unsigned"


# corFnc = ifelse(corType=="pearson", cor, bicor)
# 对二元变量，如样本性状信息计算相关性时，
# 或基因表达严重依赖于疾病状态时，需设置下面参数
# maxPOutliers = ifelse(corType=="pearson",1,0.05)

# 关联样品性状的二元变量时，设置
# robustY = ifelse(corType=="pearson",T,F)


dataExpr<- exprSet_test
dataExpr<- exprSet
dim(dataExpr)
head(dataExpr)[,1:8]

######################## step2数据筛选
## 筛选中位绝对偏差前75%的基因，至少MAD大于0.01
## 筛选后会降低运算量，也会失去部分信息
## 也可不做筛选，使MAD大于0即可
# m.mad <- apply(dataExpr,1,mad)
# dataExprVar <- dataExpr[which(m.mad > max(quantile(m.mad, probs=seq(0, 1, 0.25))[2],0.01)),]
# 
# ## 转换为样品在行，基因在列的矩阵
# dataExpr <- as.data.frame(t(dataExprVar))
# 


## 因为WGCNA针对的是基因进行聚类，而一般我们的聚类是针对样本用hclust即可，所以这个时候需要转置。
WGCNA_matrix = t(dataExpr[order(apply(dataExpr,1,mad), decreasing = T)[1:5000],])
datExpr0 <- WGCNA_matrix  ## top 5000 mad genes
dataExpr <- datExpr0 


## 检测缺失值
gsg = goodSamplesGenes(dataExpr, verbose = 3)

if (!gsg$allOK){
  # Optionally, print the gene and sample names that were removed:
  if (sum(!gsg$goodGenes)>0) 
    printFlush(paste("Removing genes:", 
                     paste(names(dataExpr)[!gsg$goodGenes], collapse = ",")));
  if (sum(!gsg$goodSamples)>0) 
    printFlush(paste("Removing samples:", 
                     paste(rownames(dataExpr)[!gsg$goodSamples], collapse = ",")));
  # Remove the offending genes and samples from the data:
  dataExpr = dataExpr[gsg$goodSamples, gsg$goodGenes]
}
nGenes = ncol(dataExpr)
nSamples = nrow(dataExpr)

dim(dataExpr)

########################### step2 软阈值筛选

## 查看是否有离群样品
sampleTree = hclust(dist(dataExpr), method = "average")
plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="")



# 软阈值的筛选原则是使构建的网络更符合无标度网络特征。

powers = c(c(1:10), seq(from = 12, to=30, by=2))
sft = pickSoftThreshold(dataExpr, powerVector=powers, 
                        networkType=type, verbose=5)

par(mfrow = c(1,2))
cex1 = 0.9
# 横轴是Soft threshold (power)，纵轴是无标度网络的评估参数，数值越高，
# 网络越符合无标度特征 (non-scale)
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")

# 筛选标准。R-square=0.85
abline(h=0.85,col="red")

# Soft threshold与平均连通性
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")
power = sft$powerEstimate
power


# ### 经验power (无满足条件的power时选用)
# # 无向网络在power小于15或有向网络power小于30内，没有一个power值可以使
# # 无标度网络图谱结构R^2达到0.8，平均连接度较高如在100以上，可能是由于
# # 部分样品与其他样品差别太大。这可能由批次效应、样品异质性或实验条件对
# # 表达影响太大等造成。可以通过绘制样品聚类查看分组信息和有无异常样品。
# # 如果这确实是由有意义的生物变化引起的，也可以使用下面的经验power值。
# if (is.na(power)){
#   power = ifelse(nSamples<20, ifelse(type == "unsigned", 9, 18),
#                  ifelse(nSamples<30, ifelse(type == "unsigned", 8, 16),
#                         ifelse(nSamples<40, ifelse(type == "unsigned", 7, 14),
#                                ifelse(type == "unsigned", 6, 12))       
#                  )
#   )
# }





############## step3 网络构建
##一步法网络构建：One-step network construction and module detection##
# power: 上一步计算的软阈值
# maxBlockSize: 计算机能处理的最大模块的基因数量 (默认5000)；
#  4G内存电脑可处理8000-10000个，16G内存电脑可以处理2万个，32G内存电脑可
#  以处理3万个
#  计算资源允许的情况下最好放在一个block里面。
# corType: pearson or bicor
# numericLabels: 返回数字而不是颜色作为模块的名字，后面可以再转换为颜色
# saveTOMs：最耗费时间的计算，存储起来，供后续使用
# mergeCutHeight: 合并模块的阈值，越大模块越少


###### 20的内存都不够用，所以maxBlockSize只能选取前5000gene
maxPOutliers=0
### 软阈值
power=sft$powerEstimate
# 相关性计算
# 官方推荐 biweight mid-correlation & bicor
# corType: pearson or bicor
# 为与原文档一致，故未修改
corType = "pearson"
exprMat <- "WGCNA/LiverFemaleClean.txt"  
##### 这个很关键,需要创建WGCNA的文件夹，然后创建LiverFemaleClean的文本
#否则报错 Error in gzfile(file, "wb") : cannot open the connection


net = blockwiseModules(dataExpr, power = power,
                       # maxBlockSize = nGenes,
                       maxBlockSize = 6000,
                       TOMType = type, minModuleSize = 30,
                       reassignThreshold = 0, mergeCutHeight = 0.25,
                       numericLabels = TRUE, pamRespectsDendro = FALSE,
                       saveTOMs=TRUE, corType = corType,
                       maxPOutliers=maxPOutliers, loadTOMs=TRUE,
                       saveTOMFileBase = paste0(exprMat, ".tom"),  ##### 这个很关键,需要创建WGCNA的文件夹，然后创建LiverFemaleClean的文本
                       verbose = 3)
### 报错 Error in (new("standardGeneric", .Data = function (x, y = NULL, use = "everything",  :    unused arguments (weights.x = NULL, weights.y = NULL, cosine = FALSE)
## 解决方案：I think I tackle the problem. There is a conflict between the WGCNA and the other packages. the other package have a function the same as the other in WGCNA in the run r studio. when I library no packages other than WGCNA, the program runs well and get the module. Thanks for professor Kevin Blighe's advice.
## 也就是不要加载别的程序，cor这个命令可能和别的包相互冲突https://www.biostars.org/p/305714/

# jimi版http://www.bio-info-trainee.com/2535.html
# net = blockwiseModules(
#   datExpr,
#   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
# )

# 根据模块中基因数目的多少，降序排列，依次编号为 `1-最大模块数`。
# **0 (grey)**表示**未**分入任何模块的基因。 
table(net$colors)

#### step4  层级聚类树展示各个模块
## 灰色的为**未分类**到模块的基因。
# Convert labels to colors for plotting
moduleLabels = net$colors
moduleColors = labels2colors(moduleLabels)
# Plot the dendrogram and the module colors underneath
# 如果对结果不满意，还可以recutBlockwiseTrees，节省计算时间
plotDendroAndColors(net$dendrograms[[1]], moduleColors[net$blockGenes[[1]]],
                    "Module colors",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)

#### step5 绘制模块之间相关性图

# module eigengene, 可以绘制线图，作为每个模块的基因表达趋势的展示
MEs = net$MEs

### 不需要重新计算，改下列名字就好
### 官方教程是重新计算的，起始可以不用这么麻烦
MEs_col = MEs
library(stringr)
colnames(MEs_col) = paste0("ME", labels2colors(
  as.numeric(str_replace_all(colnames(MEs),"ME",""))))
MEs_col = orderMEs(MEs_col)

# 根据基因间表达量进行聚类所得到的各模块间的相关性图
# marDendro/marHeatmap 设置下、左、上、右的边距
plotEigengeneNetworks(MEs_col, "Eigengene adjacency heatmap", 
                      marDendro = c(3,3,2,4),
                      marHeatmap = c(3,4,2,2), plotDendrograms = T, 
                      xLabelsAngle = 90)

# ## 如果有表型数据，也可以跟ME数据放一起，一起出图
# sampleName = rownames(dataExpr)
# pdata[1:4,1:4]
# # traitData<- data.frame(sex=pdata$`Sex:ch1`,status = pdata$`status:ch1`)
# traitData<- data.frame(sex=as.numeric(factor(pdata$`Sex:ch1`)),status = as.numeric(factor(pdata$`status:ch1`)))
# row.names(traitData)<- row.names(pdata)
# traitData[1:2,1:2]
# match(sampleName, rownames(traitData))
# traitData = traitData[match(sampleName, rownames(traitData)), ]
# MEs_colpheno = orderMEs(cbind(MEs_col, traitData))
# plotEigengeneNetworks(MEs_colpheno, "Eigengene adjacency heatmap", 
#                       marDendro = c(3,3,2,4),
#                       marHeatmap = c(3,4,2,2), plotDendrograms = T, 
#                       xLabelsAngle = 90)

#### step6 可视化基因网络 (TOM plot)

# 如果采用分步计算，或设置的blocksize>=总基因数，直接load计算好的TOM结果
# 否则需要再计算一遍，比较耗费时间
# TOM = TOMsimilarityFromExpr(dataExpr, power=power, corType=corType, networkType=type)
load(net$TOMFiles[1], verbose=T)

## Loading objects:
##   TOM

TOM <- as.matrix(TOM)
TOM[1:4,1:4]
dissTOM = 1-TOM
# Transform dissTOM with a power to make moderately strong 
# connections more visible in the heatmap
plotTOM = dissTOM^7
# Set diagonal to NA for a nicer plot
diag(plotTOM) = NA
# Call the plot function

# 这一部分特别耗时，行列同时做层级聚类
moduleColors = labels2colors(moduleLabels)


TOMplot(plotTOM, net$dendrograms, moduleColors, 
        main = "Network heatmap plot, all genes")
####### 又报错了可能是因为样本量和基因数量太大了，当我减少基因数量和样本量时，没有再次报错
# Error in TOMplot(plotTOM, net$dendrograms, moduleColors, main = "Network heatmap plot, all genes") : 
#   ERROR: number of color labels does not equal number of nodes in dissim.
# nNodes != dim(dissim)[[1]] 
#### 计算了1个小时，最好的选择是，随机挑选400个基因进行计算 
# #然后随机选取部分基因作图
# nSelect = 400
# # For reproducibility, we set the random seed
# set.seed(10);
# 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)
# # Taking the dissimilarity to a power, say 10, makes the plot more informative by effectively changing
# # the color palette; setting the diagonal to NA also improves the clarity of the plot
# plotDiss = selectTOM^7;
# diag(plotDiss) = NA;
# TOMplot(plotDiss, selectTree, selectColors, main = "Network heatmap plot, selected genes")


# step7  导出网络用于Cytoscape
probes = colnames(dataExpr)
dimnames(TOM) <- list(probes, probes)

# Export the network into edge and node list files Cytoscape can read
# threshold 默认为0.5, 可以根据自己的需要调整，也可以都导出后在
# cytoscape中再调整
cyt = exportNetworkToCytoscape(TOM,
                               edgeFile = paste(exprMat, ".edges.txt", sep=""),
                               nodeFile = paste(exprMat, ".nodes.txt", sep=""),
                               weighted = TRUE, threshold = 0,
                               nodeNames = probes, nodeAttr = moduleColors)


### step8 关联表型数据

# trait <- "WGCNA/TraitsClean.txt"
# # 读入表型数据，不是必须的
# if(trait != "") {
#   traitData <- read.table(file=trait, sep='\t', header=T, row.names=1,
#                           check.names=FALSE, comment='',quote="")
#   sampleName = rownames(dataExpr)
#   traitData = traitData[match(sampleName, rownames(traitData)), ]
# }

# sampleName = rownames(dataExpr)
# pdata[1:4,1:4]
# # traitData<- data.frame(sex=pdata$`Sex:ch1`,status = pdata$`status:ch1`)
# traitData<- data.frame(sex=as.numeric(factor(pdata$`Sex:ch1`)),status = as.numeric(factor(pdata$`status:ch1`)))
# row.names(traitData)<- row.names(pdata)
# traitData[1:2,1:2]
# match(sampleName, rownames(traitData))
# traitData = traitData[match(sampleName, rownames(traitData)), ]

### 模块与表型数据关联
if (corType=="pearson") {
  modTraitCor = cor(MEs_col, traitData, use = "p")
  modTraitP = corPvalueStudent(modTraitCor, nSamples)
} else {
  modTraitCorP = bicorAndPvalue(MEs_col, traitData, robustY=robustY)
  modTraitCor = modTraitCorP$bicor
  modTraitP   = modTraitCorP$p
}

## Warning in bicor(x, y, use = use, ...): bicor: zero MAD in variable 'y'.
## Pearson correlation was used for individual columns with zero (or missing)
## MAD.

# signif表示保留几位小数
textMatrix = paste(signif(modTraitCor, 2), "\n(", signif(modTraitP, 1), ")", sep = "")
dim(textMatrix) = dim(modTraitCor)
labeledHeatmap(Matrix = modTraitCor, xLabels = colnames(traitData), 
               yLabels = colnames(MEs_col), 
               cex.lab = 0.5, 
               ySymbols = colnames(MEs_col), colorLabels = FALSE, 
               colors = blueWhiteRed(50), 
               textMatrix = textMatrix, setStdMargins = FALSE, 
               cex.text = 0.5, zlim = c(-1,1),
               main = paste("Module-trait relationships"))

#### step9 模块内基因与表型数据关联, 

## 从上图可以看到MEmagenta与Insulin_ug_l相关

## 模块内基因与表型数据关联

# 性状跟模块虽然求出了相关性，可以挑选最相关的那些模块来分析，
# 但是模块本身仍然包含非常多的基因，还需进一步的寻找最重要的基因。
# 所有的模块都可以跟基因算出相关系数，所有的连续型性状也可以跟基因的表达
# 值算出相关系数。
# 如果跟性状显著相关基因也跟某个模块显著相关，那么这些基因可能就非常重要
# 。

### 计算模块与基因的相关性矩阵

if (corType=="pearson") {
  geneModuleMembership = as.data.frame(cor(dataExpr, MEs_col, use = "p"))
  MMPvalue = as.data.frame(corPvalueStudent(
    as.matrix(geneModuleMembership), nSamples))
} else {
  geneModuleMembershipA = bicorAndPvalue(dataExpr, MEs_col, robustY=robustY)
  geneModuleMembership = geneModuleMembershipA$bicor
  MMPvalue   = geneModuleMembershipA$p
}


# 计算性状与基因的相关性矩阵

## 只有连续型性状才能进行计算，如果是离散变量，在构建样品表时就转为0-1矩阵。

if (corType=="pearson") {
  geneTraitCor = as.data.frame(cor(dataExpr, traitData, use = "p"))
  geneTraitP = as.data.frame(corPvalueStudent(
    as.matrix(geneTraitCor), nSamples))
} else {
  geneTraitCorA = bicorAndPvalue(dataExpr, traitData, robustY=robustY)
  geneTraitCor = as.data.frame(geneTraitCorA$bicor)
  geneTraitP   = as.data.frame(geneTraitCorA$p)
}

## Warning in bicor(x, y, use = use, ...): bicor: zero MAD in variable 'y'.
## Pearson correlation was used for individual columns with zero (or missing)
## MAD.

# 最后把两个相关性矩阵联合起来,指定感兴趣模块进行分析
module = "magenta"
pheno = "sex"
modNames = substring(colnames(MEs_col), 3)
# 获取关注的列
module_column = match(module, modNames)
pheno_column = match(pheno,colnames(traitData))
# 获取模块内的基因
moduleGenes = moduleColors == module

sizeGrWindow(7, 7)
par(mfrow = c(1,1))
# 与性状高度相关的基因，也是与性状相关的模型的关键基因
verboseScatterplot(abs(geneModuleMembership[moduleGenes, module_column]),
                   abs(geneTraitCor[moduleGenes, pheno_column]),
                   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)



#### step10 分步法展示每一步都做了什么
## 计算邻接矩阵
adjacency = adjacency(dataExpr, power = power)

### 把邻接矩阵转换为拓扑重叠矩阵，以降低噪音和假相关，获得距离矩阵。
TOM = TOMsimilarity(adjacency)
dissTOM = 1-TOM

### 层级聚类计算基因之间的距离树 
geneTree = hclust(as.dist(dissTOM), method = "average")

### 模块合并
# We like large modules, so we set the minimum module size relatively high:
minModuleSize = 30
# Module identification using dynamic tree cut:
dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM,
                            deepSplit = 2, pamRespectsDendro = FALSE,
                            minClusterSize = minModuleSize)
# Convert numeric lables into colors
dynamicColors = labels2colors(dynamicMods)

### 通过计算模块的代表性模式和模块之间的定量相似性评估，合并表达图谱相似的模块
MEList = moduleEigengenes(dataExpr, colors = dynamicColors)
MEs = MEList$eigengenes
# Calculate dissimilarity of module eigengenes
MEDiss = 1-cor(MEs)
# Cluster module eigengenes
METree = hclust(as.dist(MEDiss), method = "average")
MEDissThres = 0.25

# Call an automatic merging function
merge = mergeCloseModules(dataExpr, dynamicColors, cutHeight = MEDissThres, verbose = 3)
# The merged module colors
mergedColors = merge$colors;
# Eigengenes of the new merged

## 分步法完结