使用pbmc数据集做演示，使用矩阵是@data矩阵，基因是2000高变基因，细胞是所有细胞。（pbmc的数据不适合拿来跑WGCNA，这里只做为演示。）
1. 导入数据
library(WGCNA)
pbmc <- readRDS("pbmc.rds")
datadf <- as.matrix(pbmc@assays$RNA@data )
dim(datadf)
df <- datadf[intersect(Seurat::VariableFeatures(pbmc),rownames(datadf)),]
dim(df)
df[1:4,1:4]
class(df)
# [1] "matrix" "array"
# 转换为样品在行，基因在列的矩阵
dataExpr <- as.data.frame(t(df))
dim(dataExpr)
head(dataExpr)[,1:8]
nGenes = ncol(dataExpr)
nSamples = nrow(dataExpr)
WGCNA基本参数设置
type = "unsigned" # 官方推荐 "signed" 或 "signed hybrid"
corType = "pearson" # 相关性计算 官方推荐 biweight mid-correlation & bicor corType: pearson or bicor
corFnc = ifelse(corType=="pearson", cor, bicor)
corFnc
maxPOutliers = ifelse(corType=="pearson",1,0.05) # 对二元变量，如样本性状信息计算相关性时， # 或基因表达严重依赖于疾病状态时，需设置下面参数
# 关联样品性状的二元变量时，设置
robustY = ifelse(corType=="pearson",T,F)
2. 检测缺失值和离群样本
#检测缺失值
gsg = goodSamplesGenes(dataExpr, verbose = 3)
gsg$allOK
gsg$goodSamples
# 查看是否有离群样品
sampleTree = hclust(dist(dataExpr), method = "average")
plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="")
3. 挑选软阈值
powers = c(c(1:10), seq(from = 12, to=30, by=2))
sft = pickSoftThreshold(dataExpr, powerVector=powers,
networkType="signed", 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")
参数beta取值默认是1到30，上述图形的横轴均代表权重参数β，左图纵轴代表对应的网络中log(k)与log(p(k))相关系数的平方。相关系数的平方越高，说明该网络越逼近无网路尺度的分布。右图的纵轴代表对应的基因模块中所有基因邻接函数的均值。最佳的beta值就是sft$powerEstimate，已经被保存到变量了，不需要知道具体是什么，后面的代码都用这个即可，在本例子里面是6。
其他方法挑选power值
# 自动挑选软阈值
power = sft$powerEstimate
softPower = power
softPower
# [1] 6
# 经验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))
)
)
}
power
# [1] 6
4. 构建共表达网络
#一步法网络构建：One-step network construction and module detection##
cor <- WGCNA::cor
net = blockwiseModules(dataExpr, power = power, maxBlockSize = nGenes,#nGenes
TOMType = "unsigned", minModuleSize = 10,
reassignThreshold = 0, mergeCutHeight = 0.25,
numericLabels = TRUE, pamRespectsDendro = FALSE,
saveTOMs=TRUE, corType = corType,
maxPOutliers=maxPOutliers, loadTOMs=TRUE,
saveTOMFileBase = paste0("dataExpr", ".tom"),
verbose = 3)
# power: 上一步计算的软阈值
# maxBlockSize: 计算机能处理的最大模块的基因数量 (默认5000)；
# 4G内存电脑可处理8000-10000个，16G内存电脑可以处理2万个，32G内存电脑可
# 以处理3万个
# 计算资源允许的情况下最好放在一个block里面。
# corType: pearson or bicor
# numericLabels: 返回数字而不是颜色作为模块的名字，后面可以再转换为颜色
# saveTOMs：最耗费时间的计算，存储起来，供后续使用
# mergeCutHeight: 合并模块的阈值，越大模块越少
# type = unsigned
计算过程
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 not use multithreading.
Fraction of slow calculations: 0.000000
..connectivity..
..matrix multiplication (system BLAS)..
..normalization..
..done.
..saving TOM for block 1 into file dataExpr.tom-block.1.RData
....clustering..
....detecting modules..
....calculating module eigengenes..
....checking kME in modules..
..removing 1 genes from module 2 because their KME is too low.
..merging modules that are too close..
mergeCloseModules: Merging modules whose distance is less than 0.25
Calculating new MEs...
查看一下返回的net结果
table(net$colors)
# 0 1 2 3 4
# 1855 59 52 17 17
net$unmergedColors
head(net$MEs)
# ME4 ME2 ME1 ME3 ME0
# AAACATACAACCAC-1 0.004299074 -0.0024229364 -0.012777424 -0.012259802 -0.011847288
# AAACATTGAGCTAC-1 -0.007987395 -0.0024229364 -0.007647816 0.036779648 0.003124504
# AAACATTGATCAGC-1 -0.005377686 -0.0024229364 -0.010071728 -0.016748944 -0.001960111
# AAACCGTGCTTCCG-1 -0.005498834 -0.0002653717 0.037286420 0.018902740 0.027428832
# AAACCGTGTATGCG-1 0.053098755 -0.0024229364 -0.009787707 -0.007069281 -0.002900237
# AAACGCACTGGTAC-1 -0.010385161 -0.0024229364 -0.008945768 -0.005881210 -0.008029291
net$goodSamples
net$goodGenes
net$TOMFiles
# [1] "dataExpr.tom-block.1.RData" "dataExpr.tom-block.2.RData" "dataExpr.tom-block.3.RData"
net$blockGenes
net$blocks
net$MEsOK
net$MEs是所有细胞*5个模块的矩阵（可以类比PCA的PC来理解）。⚠️这个矩阵是非常重要的，后面做随机森林也是用的这个矩阵。
这里用不同的颜色来代表那些所有的模块，其中灰色默认是无法归类于任何模块的那些基因，如果灰色模块里面的基因太多，那么前期对表达矩阵挑选基因的步骤可能就不太合适。
## 灰色的为**未分类**到模块的基因。
# Convert labels to colors for plotting
moduleLabels = net$colors
moduleColors = labels2colors(moduleLabels)
moduleColors
# 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)
哇哦，基本全是灰的
灰色模块一般如果基因不多的话，在后续分析的时候去掉就可以了。但是这个里面太多了，应该是前面基因的选择有点问题，还有就是用所有细胞来做的WGCNA，可能矩阵还是太稀疏了。
5. 网络可视化 (TOM plot)
# module eigengene, 可以绘制线图，作为每个模块的基因表达趋势的展示
MEs = net$MEs
### 不需要重新计算，改下列名字就好
### 官方教程是重新计算的，起始可以不用这么麻烦
MEs_col = MEs
colnames(MEs_col) = paste0("ME", labels2colors(
as.numeric(str_replace_all(colnames(MEs),"ME",""))))
MEs_col = orderMEs(MEs_col)
# 根据基因间表达量进行聚类所得到的各模块间的相关性图
# marDendro/marHeatmap 设置下、左、上、右的边距
head(MEs_col)
?plotEigengeneNetworks
plotEigengeneNetworks(MEs, "Eigengene adjacency heatmap",
marDendro = c(3,3,2,4),
marHeatmap = c(3,4,2,2),
plotDendrograms = T,
xLabelsAngle = 90)
6. 模块-表型数据关联
datTraits=pbmc@meta.data
MEs0 = moduleEigengenes(dataExpr, moduleColors)$eigengenes
MEs = orderMEs(MEs0)
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(mar = c(6, 8.5, 3, 3));
# 用热图的形式展示相关系数
labeledHeatmap(Matrix = moduleTraitCor,
xLabels = names(datTraits),
yLabels = names(MEs),
ySymbols = names(MEs),
colorLabels = FALSE,
colors = greenWhiteRed(50),
textMatrix = textMatrix,
setStdMargins = FALSE,
cex.text = 0.5,
zlim = c(-1,1),
main = paste("Module-trait relationships"))
灰色太多了，换个数据可能好点
7. 随机森林筛选重要的模块
identical(rownames(MEs_col),rownames(datTraits))
exp <- cbind(MEs_col,datTraits)
View(exp)
pbmc.rf <- randomForest(cell_type ~ ., data=exp, importance=TRUE,
proximity=TRUE)
library(randomForest)
pbmc.rf <- randomForest(cell_type ~ ., data=exp, importance=TRUE,
proximity=TRUE)
print(pbmc.rf)
pbmc.rf$importance
varImpPlot(pbmc.rf, main = "Top 30 - variable importance")
因为我把metadata都放进来了，所以筛出来最有意义的是cluster分群结果（毕竟细胞类型注释就是基于cluster做的）。只看几个模块的话，最重要的是模块3。
思路
1. 需要注意的是，整合样本做WGCNA分析分析得到的模块主要和细胞类型有关联，因为此时的cluster代表不同的细胞类型，此时做WGCNA感觉作用不大，因为不同的细胞类型本身就应该有独特的一套基因列表，但是我们在做同一细胞类型再分群的时候，这个分析的用处就会很大，因为得到的模块与同一细胞类型不同的subcluster相关联，得到的就是一些在各个亚群高度协调的模块基因，这个时候把得到的模块基因全部做下游分析，意义非常大，从另外一个侧面表现了细胞类型内部的异质性，如果是多样本整合的关系（比如正常和疾病），同一细胞类型的不同样本关联得到的gene模块列表，就可以得到处理组和对照组关联度高的模块列表，新的角度来看待疾病的发生。
2. 此外，我们在做分析的时候，同样是多个样本（正常和疾病两组），同一细胞类型，我们先做正常组（或者疾病组）的WGCNA，得到的基因模块反映到疾病组（正常组），这个时候就会看出来明显的差异，本应该高度关联的模块可能就消失了，这个高度关联的模块的功能也随之产生了变化；还有就是同一细胞类型，两组样本先分别做WGCNA，得到的基因模块进行平行对比，那么得到的就是对疾病发生新的认识。
