首先，WGCNA并不适合在2.1的差异基因筛选后的基因集上做，原因在知乎上有大佬们进行过讨论：
@云生信学生物信息：
强烈不建议采用差异表达的基因进行WGCNA分析， 因为差异表达的基因，就是的样本失去了多样性。
@Zheng博士：
WGCNA旨在识别共表达的基因模块，而不是单个基因，分析结果应在此背景下进行解释。使用WGCNA来鉴定DEG可能没有考虑到模块内基因之间的潜在相互作用，并且可能会错过重要的生物学见解
@胡喵喵
WGCNA求的是基因共表达网络，用筛选后的差异表达基因分析，网络节点变少，会降低稳定性。
一些其他的参考链接：
TCGA数据挖掘实战（二）——WGCNA①：
不建议直接使用差异基因进行分析，因为差异基因本身就会形成单个或几个高度相关的模块，与WGCNA设计冲突。
小鬼的WGCNA分析详解(三)-挖掘与表型(体重)相关的基因 - 简书 (jianshu.com)
本文作者推荐了一个参考文献：Integrating Genetic and Network Analysis to Characterize Genes Related to Mouse Weight；
“官方给的是不建议使用差异表达基因做：”
“文章中基因的筛选策略为： 对于模块检测，我们将分析限制在 3,600 个连接最紧密的基因上，因为我们的模块构建方法和可视化工具目前无法处理更大的数据集。根据定义，模块基因与其模块的基因高度连接（即模块基因往往具有相对较高的连接性）。因此，出于模块检测的目的，将分析限制在连接最紧密的基因上不应导致重大信息丢失。由于我们分析中的网络节点对应于基因而不是探针集，因此我们消除了同一基因具有相似表达模式的多个探针。 ”【简书作者引用的翻译】
其次，讲述WGCNA的原理：
A1：WGCNA分析是什么方法
WGCNA（Weighted Gene Co-expression Network Analysis）是一种系统生物学方法，用于分析基因表达数据中基因间的共表达模式，构建基因共表达网络，并识别模块化的基因集合。它着重于探索基因之间的关联性，而不是单个基因的差异性。主要步骤包括构建基因共表达网络、识别模块、发现模块与临床特征相关性等。
A2：针对什么问题，得到什么结论，运用在哪个方面？
问题：WGCNA主要针对基因表达数据中的模式发现和相关性探索，能够揭示不同基因之间的共表达关系，从而找出共同调控的基因模块。结论：WGCNA可以识别具有相似表达模式的基因集合，生成基因模块并发现与特定性状（如疾病状态、临床特征等）相关的模块。应用方面：WGCNA广泛应用于生物医学研究，包括疾病机制的探索、生物标志物的鉴定、疾病分类、疾病预后等。在乳腺癌研究中，WGCNA可用于发现与乳腺癌相关的基因模块，以及这些模块与乳腺癌生存期、病理特征等的关联。
A3：WGCNA分析在建立预后模型中的作用？
WGCNA在乳腺癌预后模型建立中的作用：WGCNA用于乳腺癌预后模型建立的初步步骤，以识别具有共表达模式的基因模块，并探索这些模块与临床特征的关系。它能够帮助筛选出与乳腺癌生存期相关的基因集合。
代码：参考：
TCGA数据挖掘实战（二）——WGCNA①
step1：1.样本聚类图和性状热图
知识点：样本聚类图和形状热图
注意，这个时候做样本聚类只是为了剔除离群值。
# 剔除离群样本 sampleTree <- hclust(dist(datExpr0), method = "average") # 这段代码中使用了hclust函数，它是R语言中用于进行层次聚类分析的工具之一。 # dist函数计算了基因表达数据集datExpr0的距离矩阵，然后hclust使用所计算的距离矩阵进行层次聚类分析，方法选用了"average"（平均连接）的方式。 # 这个步骤的目的是将样本进行聚类，形成一个层次聚类树，以便于发现样本之间的相似性和差异性。
关键点：
使用的是第一步里的exprSet先取log（其实前面还存了exprSet_log可以直接用）
对样本进行了质量控制，剔除了可能是离群值的样本，然后利用聚类分析（hierarchical clustering）对样本进行聚类，以便观察样本之间的关系，并剔除聚类结果中的离群样本。
基础筛选：使用了 apply 函数对 exprData 中的每一行进行计算，计算每行数据的MAD（中位数绝对偏差）。根据MAD对 exprData 进行筛选，保留变异性较大的前5000行数据。用了goodSamplesGenes()函数继续做筛选：# goodSamplesGenes() 函数通常用于检查数据中是否存在质量较差的样本或基因，并对它们进行过滤或移除，以确保后续分析的准确性和可靠性。（在这里，该函数用于评估和识别具有较多缺失值或低表达水平的样本和基因。）提出离群样本：使用层次聚类的方法，将明显离群的树切除（这里设定的高度是235）
################### 1.样本聚类图和性状热图 ############################
setwd("D:\\Datamining3")
options(stringsAsFactors = F)
library(WGCNA)
library(readr)
counts <- read.table("BRCA_expr.txt", header = T, check.names = F)
# 读取counts矩阵
exprData <- log2(edgeR::cpm(counts+1))
# 对其进行log2(cpm+1)操作之后可以直接用于WGCNA分析。
# 上述操作对数据进行了对数转换，使得数据更接近正态分布，且避免了零值所带来的问题。
# 其中CPM代表Counts Per Million，是一种常见的标准化方法，用于调整不同样本之间的大小差异。
head(exprData)
# 筛选数据
exprMad <- apply(exprData, 1, mad)
# 这段代码使用了 apply 函数对 exprData 中的每一行进行计算，计算每行数据的MAD（中位数绝对偏差）。
datExpr0 <- as.data.frame(t(exprData[order(exprMad, decreasing = T)[1:5000],]))
# 根据MAD对 exprData 进行筛选，保留变异性较大的前5000行数据。
# check for genes and samples with too many missing values
gsg <- goodSamplesGenes(datExpr0, verbose = 3);
# goodSamplesGenes() 函数通常用于检查数据中是否存在质量较差的样本或基因，并对它们进行过滤或移除，以确保后续分析的准确性和可靠性。
# 在这里，该函数用于评估和识别具有较多缺失值或低表达水平的样本和基因。
gsg$allOK
# TRUE
# 剔除离群样本
sampleTree <- hclust(dist(datExpr0), method = "average")
# 这段代码中使用了hclust函数，它是R语言中用于进行层次聚类分析的工具之一。
# dist函数计算了基因表达数据集datExpr0的距离矩阵，然后hclust使用所计算的距离矩阵进行层次聚类分析，方法选用了"average"（平均连接）的方式。
# 这个步骤的目的是将样本进行聚类，形成一个层次聚类树，以便于发现样本之间的相似性和差异性。
par(cex = 0.5)
par(mar = c(0,6,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 = 235, col = "Red") #
# Determine cluster under the line
clust = cutreeStatic(sampleTree, cutHeight = 235, minSize = 10)
# cutreeStatic`函数用于根据先前构建的样本聚类树`sampleTree`以某个高度进行切割
# 并且设定每个簇（cluster）最少包含的样本数为10（`minSize = 10`）。
table(clust)# 列出切割后各个簇中样本的数量，感觉比较合适
outSamples <- datExpr0[(clust!=1),]# (clust == 0)`选取被划分到第一个簇的样本，将它们从原始数据集`datExpr0`中提取出来，并赋值给`outSamples`。
rownames(outSamples)
# 最后，`rownames(outSamples)`返回`outSamples`中行的名称，即这些样本的标识。
# clust 1 contains the samples we want to keep.
keepSamples <- (clust==1)
datExpr <- datExpr0[keepSamples, ]
# 根据指定的原则keepsamples来挑选数据
sampleTree2 <- hclust(dist(datExpr), method = "average");
par(cex = 0.5);
par(mar = c(0,6,2,0))
plot(sampleTree2, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5,
cex.axis = 1.5, cex.main = 2)
加入临床特征，就开始做WGCNA分析了：
关键点：
临床特征有很多很多，这里能作为输入的必须为数值型特征。这里选了4个：
"submitter_id.samples","age_at_initial_pathologic_diagnosis", "number_of_lymphnodes_positive_by_he", "days_to_death.demographic","initial_weight.samples"
- **submitter_id.samples：** 这是样本的提交者标识符，用于唯一标识每个样本或数据提交者。 "age_at_initial_pathologic_diagnosis"（初次病理诊断年龄）是指患者首次被进行病理学诊断时的年龄。在癌症研究中，这个特征通常表示患者首次被诊断为患有某种癌症时的年龄。这个特征可以用来研究年龄与癌症发生的关系，比如探索癌症的发病年龄分布、不同年龄段发病风险的差异等。对于分析癌症的研究，这个特征可以作为一个重要的预测因子或者研究因素。 "number_of_lymphnodes_positive_by_he" 表示淋巴结中受到肿瘤影响的阳性淋巴结数量。 initial_weight.samples" 可能是描述研究样本在开始时的体重或重量。【患者整体健康状况】在医学研究中，初始体重信息可能被用来分析患者的健康状况、疾病发展和治疗反应等方面。 **days_to_death.demographic：** 这个字段表示从特定时间点（如诊断日期或其他特定事件）到患者死亡的天数。它提供了患者生存时间或存活期限的信息。 此处的临床特征选择参考：基于R对TCGA患者临床数据整理和Cox分析 - 简书 一般常规项包括：submitter_id，vital_status，tumor_stage，pathologic_M，pathologic_N，pathologic_T，age_at_initial_pathologic_diagnosis，days_to_death，days_to_last_follow_up，race，gender但具体选择哪些临床特征还需根据研究的具体问题和数据的可用性来决定。一般常规项包括：submitter_id，vital_status，tumor_stage，pathologic_M，pathologic_N，pathologic_T，age_at_initial_pathologic_diagnosis，days_to_death，days_to_last_follow_up，race，gender但具体选择哪些临床特征还需根据研究的具体问题和数据的可用性来决定。画个图：样本聚类树状图(dendrogram)和性状热图(traitheatmap)，展示了样本聚类和临床特征之间的关系。
# 加入临床特征（性状）
phenotype <- read_tsv("TCGA-BRCA.GDC_phenotype.tsv.gz")
dim(phenotype)
names(phenotype)
colnames(phenotype)
final_expr.GBM[1:6,36:39]
##看看样本里的分类
sample <- ifelse(substring(phenotype$submitter_id.samples,14,15)<10, "tumor", "normal")
sample <- factor(sample, levels = c("tumor", "normal"))
group <- data.frame(barcode = phenotype$submitter_id.samples, sample)
table(group$sample)
numerics <- sapply(phenotype, class)
# 检查结果并存储数值型变量
numeric_variables <- names(numerics[numerics == "numeric"])
numeric_variables
# 提取 submitter_id.samples 列
submitter_id <- phenotype[["submitter_id.samples"]]
# 提取 numeric_variables 中指定的列
numeric_cols <- phenotype[, numeric_variables]
# 将提取的列添加到 allTraits 中
albind(submitter_id, numeric_cols)
allTraits
allTraits <- phenotype[,c("submitter_id.samples","age_at_initial_pathologic_diagnosis",
"number_of_lymphnodes_positive_by_he", "days_to_death.demographic","initial_weight.samples")]
###一定要是数值型
###看看自己想要的特征
allTraits <- as.data.frame(allTraits)
samples <- rownames(datExpr)
traitRows <- match(samples, allTraits$submitter_id.samples)
datTraits <- allTraits[traitRows,]
rownames(datTraits) <- datTraits$submitter_id.samples
datTraits <- datTraits[,-1]
datTraits[,3]<- as.numeric(datTraits[,3])
# 这里将样本临床特征与筛选过的样本关联起来，使样本排列顺序相同。看一下处理好的数据。
datTraits[1:5,]
sapply(phenotype, class) ##检验一下是不是数值型
# Convert traits to a color representation: white means low, red means high, grey means missing entry
traitColors <- numbers2colors(datTraits, signed = FALSE);
# Plot the sample dendrogram and the colors underneath.
plotDendroAndColors(sampleTree2, traitColors,
groupLabels = names(datTraits),
main = "Sample dendrogram and trait heatmap")
# save(datExpr, datTraits, file = "WGCNAdata/COAD-01-dataInput.RData")
plotDendroAndColors #优化一下
save(datExpr, datTraits, file = "BRCA-01-dataInput.RData")
save(datExpr, datTraits, file = "BRCA-01-dataInput_11.24.RData")
网络构建和模块检测
选择软阙值
chat老师：
下面解释一下为什么要使用软阈值以及为什么选择软阈值。
1.稀疏性： 软阈值的主要目的之一是增强网络的稀疏性。通过对共表达矩阵进行软阈值处理，将弱相关性的基因对的权重置为零，从而在网络中引入稀疏性。这样做有助于减少网络中的噪声和不相关信息，更好地捕捉真实的生物学模式。
2.强调关键关系： 软阈值化可以突出具有相对较高相关性的基因对，而将低相关性的基因对的权重降低。这有助于保留网络中那些更具生物学意义的关联，而忽略较弱的、可能是由于噪声引起的关联。
3.模块检测： WGCNA的目标之一是将基因分组为模块，即具有相似表达模式的基因的集合。软阈值化有助于明确定义模块，因为它使得只有相对较强的相关性才会形成模块，而弱相关性则被剔除。
4.网络的可解释性： 硬阈值可能导致网络过于稠密，使得难以解释。软阈值化可以使网络更加稀疏，从而减少网络的复杂性，使其更容易解释和理解。
但是发现还是一知半解，直到看到了大佬的文章：
小鬼的WGCNA图文详解(一)--软阈值 - 简书 (jianshu.com)
那么我们做WGCNA的时候为什么不用硬阈值而是选择软阈值呢？
为了说清楚这个软阈值的问题，进行一下思路的说明：
1，网络中边的理解：WGCNA分析全称为Weighted Gene Co-Expression Network Analysis，即加权的基因共表达网络分析。
网络中基因与基因之间是否连边取决于这两个基因是否发生显著共表达：nodes represent genes and nodes are connected if the corresponding genes are significantly co-expressed across appropriately chosen tissue samples
2，真实网络与随机网络的特征差异
1，网络中边的理解
计算共表达有很多方法，官网教程中给出的方法是用皮尔森相关系数(Pearson) correlation coefficient：It is standard to use the (Pearson) correlation coefficient as a co-expression measure, e.g., the absolute value of Pearson correlation is often used in a gene expression cluster analysis.
那么现在问题来了：我得到了两个基因的相关系数值之后，如何决定这两个基因在构建网络时是否连边呢？
此时你可能会想到：
1，定义一个阈值，比如有统计学意义的r>0.9，那么就连边；否则，就不连边。这个办法就是‘hard’ threshold（硬阙值），也即不加权（unweighted）。
这个办法有以下缺点：
第一：如何确定这个阈值？这个阈值比较人为，有人觉得0.9就是强相关，也有人会认为0.8也是强相关。
第二：对于真实的生物学网络，这种二元定义连边的方法是否真的适合？
2，真实网络与随机网络的特征差异
真实的生物学网络比如：
human disease network，gene co-expression networks，protein-protein interaction networks，cell-cell interaction networks，the world wide web and social interaction networks。
他们都有一个特征，那就是网络中节点的度服从幂率分布（power law） p(k)~ k−γ，即无标度网络(Scale-free networks）
而随机网络一般服从泊松(Poisson)分布
这两个网络分布的特点是：
真实网络中绝大部分点的度都很低，只有少数的点（即常说的hub节点）度很高，而随机网络中绝大部分的点的度都处于上图中的峰值处，即度相似。无标度网络节点的度分布特征使得它有一个很大的特点，那就是稳健性：随机去除网络中的一个节点，网络还能依然保持（绝大多数节点的度很低）。但是它也有脆弱的一面：那就是去除网络中的hub节点，网络就散了。但是hub节点只是网络中极少数的几个节点，被攻击的概率非常小。
基于以上，即‘hard’ threshold的缺点和应用于真实网络造成的信息丢失（参考文献中有说明）等原因，官网教程中的作者提出了一种方法即用‘soft’ thresholding 来权衡网络中的连边。
就是将相关系数的值[0,1]通过一个换算映射到[0,1]，也即加权（weighted）的思想：
β即软阙值
第二步，作者@_十三给出了两个问题：
如何选取Soft Threshold，多少的时候合适？（即上图中左右图的意思和纵坐标的理解）对于这张图的说明，官网给出的注释是这个样子的：
我做的图：
尝试用@_十三的解释来解释我的图：
【官网解释】Figure 1: Analysis of network topology for various soft-thresholding powers.
The left panel shows the scale-free fit index (y-axis) as a function of the soft-thresholding power (x-axis). 【无标度拟合指数】The right panel displays the mean connectivity(degree, y-axis) as a function of the soft-thresholding power (x-axis). 【平均连通性(度,y轴)】
【@_十三】：
1. 左图的纵坐标scale-free fit index，即signed R2，代表对应的网络中log(k)与log(p(k))相关系数的平方乘以一个方向向量，由slope决定（The sign of the scale-free model fitting index R2 is determined by minus the sign of the slope），拟合的线性方程为下图，来源于WNCGCNA包中的源代码：
相关系数的平方越高，说明该网络越逼近无标度网络的分布。相关参考文献中有大量数据证明当signed R2 大于0.85时，网络就已经符合无标度网络的分布。
因此，WGCNA包中计算SoftThreshold的函数pickSoftThreshold中RsquaredCut 默认值为 0.85，最佳的powers值保存在sft$powerEstimate。（注：因此后续R中用sft$powerEstimate挑最优，挑出来就是7）
2，右图的纵轴代表对应的网络中所有基因连接数（即节点的度）的均值。
作者：_十三 链接：https://www.jianshu.com/p/a783e06bb1a5 来源：简书 著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。
代码演示：（根据实际情况修改）TCGA数据挖掘实战（二）——WGCNA①
rm(list = ls())
options(stringsAsFactors = F)
library(WGCNA)
# 选择软阈值
# 使用函数**pickSoftThreshold()**选择适当的软阈值。
load("BRCA-01-dataInput.RData")
load("BRCA-01-dataInput_11.24.RData")
# Choose a set of soft-thresholding powers
powers <- c(c(1:10), seq(from = 12, to=20, by=2))
# Call the network topology analysis function
sft <- pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)
# Plot the results:
par(mar = c(3,6,2,1))
par(mfrow = c(1,2))
cex1 = 1
# 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", # n表示不绘图
main = paste("Scale independence"))
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
labels=powers, cex=cex1, col="red")
# best：推荐的最优软阈值，最优值为7
sft$powerEstimate
# best：推荐的最优软阈值，最优值为5（11.24版）
sft$powerEstimate
# this line corresponds to using an R^2 cut-off of h
abline(h=0.85,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")
abline(h=10,col="red")
#### 左图y轴是无标度拓扑拟合指数，右图y轴是平均连通度。
一步式网络构建和模块检测（聚类树状图）
# 注释参考：https://zhuanlan.zhihu.com/p/168262193
net <- blockwiseModules(datExpr,
power = 7, # 选前面选出来的最优阙值
TOMType = "unsigned", # 构建无尺度网络
minModuleSize = 30, # 最小模块基因数为30
reassignThreshold = 0,
mergeCutHeight = 0.25, # 模块合并阀值
numericLabels = TRUE, # 模块颜色标签
pamRespectsDendro = FALSE,
saveTOMs = TRUE, # 保存TOM矩阵
saveTOMFileBase = "BRCATOM",
verbose = 3)
#minModuleSize表示用于模块检测的最小模块尺寸。
#mergeCutHeight表示用于模块合并的树形图切割高度。这两个值越大，模块越少。
#saveTOMFileBase用来设置数据保存位置及名称。
知识点：跟着一起实验学WGCNA 最“细腻小白”的学习笔记分享（3） - 知乎 (zhihu.com)
TOM矩阵，topological overlap matrix，TOM拓扑重叠矩阵。
TOM矩阵里面的数值都在[0,1]间，是用于反应两两基因共表达的相似度，值越高，说明共表达相似度越高。即：数值越高，二者越可能在一个模块里。
结果：总共得到了6个模块
> table(net$colors)#检测出的模块数量和每个模块中包含的基因数量
0 1 2 3 4 5 6
2567 1545 402 227 139 82 38
# 查看划分的模块数和每个模块里面包含的基因个数。
# 0代表无法识别的基因数。
# 用于模块识别的分层聚类树形图存储在net$dendprograms[[1]]中。可以与模型颜色分配一起显示。
下一步，绘图：TOM图
sizeGrWindow(12, 9)
# Convert labels to colors for plotting
mergedColors <- labels2colors(net$colors)
# Plot the dendrogram and the module colors underneath
plotDendroAndColors(net$dendrograms[[1]],
mergedColors[net$blockGenes[[1]]],
"Module colors",
dendroLabels = FALSE,
hang = 0.03,
addGuide = TRUE,
guideHang = 0.05)
这个图可以分为两部分看：上半部分是基因的层次聚类树状图，下半部分是基因模块，也就是网络模块。上下对应，可以看到距离较近的基因（聚类到同一条分支）被划分到了同一模块。
每一个树枝代表一个gene，纵坐标的Heigth为聚类距离。 设定一个聚类距离阈值，就可以将gene聚成module，这里又有很多方法设定这个阈值。而官网教程成给定的是dynamicTreeCut包中的函数cutreeDynamic，并且设定了最小模块中包含的基因个数为30个基因。【本文用的是blockwiseModules函数】
官方注释：
Clustering dendrogram of genes, with dissimilarity based on topological overlap, together with assigned merged module colors and the original module colors.
### 保存参数
moduleLabels <- net$colors
moduleColors <- labels2colors(net$colors)
MEs <- net$MEs;
geneTree <- net$dendrograms[[1]];
save(MEs, moduleLabels, moduleColors, geneTree,
file = "BRCA-02-networkConstruction-auto.RData")
注意，后面还画了一个图：小鬼的WGCNA图文详解(六)-visualize a weighted network - 简书 (jianshu.com)
【官网注释】：
Visualizing the gene network using a heatmap plot. The heatmap depicts the Topological Overlap Matrix(TOM) among all genes in the analysis. Light color represents low overlap and progressively darker red color represents higher overlap. Blocks of darker colors along the diagonal are the modules. The gene dendrogram and module assignment are also shown along the left side and the top.
使用热图可视化基因网络。热图描绘了分析中所有基因之间的拓扑重叠矩阵（TOM）。浅色表示低重叠，逐渐变深的红色表示较高重叠。沿对角线的深色块是模块。基因树状图和模块分配也显示在左侧和顶部。
两个部分：1.聚类树 2. 热图
1. 聚类树（前面已经出现过了）
2. 热图：
热图的对角线可以看出有一些深颜色的模块，就是相关性大的基因聚类成的模块，每一个色块都都是一个模块，与旁边的聚类树底下的模块对应。可以很明显的看出，模块内部的基因相关性很大。除此之外，其实，还可以看出模块与模块之间也是有关联而不是相互独立的。
怎么画？（用法2，法1会卡）
# Tom热图：不要用下面这种，会卡爆，后面给了快的
rm(list = ls())
options(stringsAsFactors = F)
library(WGCNA)
load("BRCA-01-dataInput.RData")
load("BRCA-02-networkConstruction-auto.RData")
nGenes <- ncol(datExpr)
nSamples <- nrow(datExpr)
# Calculate topological overlap anew: this could be done more efficiently by saving the TOM
# calculated during module detection, but let us do it again here.
dissTOM <- 1-TOMsimilarityFromExpr(datExpr, power = 6);
# 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
sizeGrWindow(9,9)
TOMplot(plotTOM, geneTree, moduleColors, main = "Network heatmap plot, all genes")
###tom热图，简单版
#由于使用所有的基因进行绘图需要比较长的时间，因此可以对基因进行抽样，随机选择部分基因进行绘图，加快绘图速度，且得到的结果相似。
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")
TOMplot((1-plotDiss), selectTree, selectColors, main ="Network heatmap plot, selected genes")
step2：基于所有基因表达谱，利用WGCNA构建共表达网络(肿瘤和正常性状)并作图，筛选出和肿瘤最相关的模块。
模块与性状的关系
既然已经得到了模块，我们就需要评估一下模块情况。所以，下一步计算模块特征向量和临床性状之间相关系数矩阵。
# 模块与性状关系图
rm(list = ls())
options(stringsAsFactors = F)
library(WGCNA)
load("BRCA-01-dataInput.RData")
load("BRCA-02-networkConstruction-auto.RData")
# Define numbers of genes and samples
nGenes <- ncol(datExpr) # 计算了数据框 datExpr 的列数，即基因的数量
nSamples <- nrow(datExpr)# 计算了数据框 datExpr 的行数，即样本的数量
# 重新计算模块特征向量（MEs），并计算这些特征向量与表型数据的相关性及相关性的 p 值。
MEs0 <- moduleEigengenes(datExpr, moduleColors)$eigengenes
MEs0
# 通过调用 moduleEigengenes 函数，使用基因表达数据框 datExpr 和模块颜色信息 moduleColors 计算模块特征向量。
# $eigengenes 部分表示从返回的结果中提取模块特征向量
MEs <- orderMEs(MEs0)
# 通过调用 orderMEs 函数，对模块特征向量进行排序。这通常是为了更好地可视化或进一步分析。
moduleTraitCor <- cor(MEs, datTraits, use = "p");
# 通过调用 cor 函数，计算模块特征向量 MEs 与表型数据 datTraits 之间的相关性矩阵，并将结果存储在变量 moduleTraitCor 中。参数 use = "p" 表示使用 Pearson 相关系数进行计算。
moduleTraitPvalue <- corPvalueStudent(moduleTraitCor, nSamples)
# 通过调用 corPvalueStudent 函数，计算相关性矩阵 moduleTraitCor 中每个元素的 p 值这一步用于评估相关性的统计显著性。
sizeGrWindow(1,1)
# 设置图形窗口的大小。
textMatrix <- paste(signif(moduleTraitCor, 2), "\n(",
signif(moduleTraitPvalue, 1), ")", sep = "")
# 这一行代码使用 paste 函数，将模块特征向量与表型数据相关性矩阵 (moduleTraitCor) 中的数值和 p 值以特定格式组合成文本，
# 并将结果存储在变量 textMatrix 中。
# signif 函数用于限制数字的小数位数：相关阵保留两位小数，p值保留一位小数
dim(textMatrix) <- dim(moduleTraitCor)
# 通过 dim 函数，将 textMatrix 的维度设置为与相关性矩阵 moduleTraitCor 相同的维度。这确保文本矩阵与相关性矩阵的行列数相匹配。
par(mar = c(6, 10, 3, 3))
# Display the correlation values within a heatmap plot
labeledHeatmap(Matrix = moduleTraitCor,
xLabels = names(datTraits),
yLabels = names(MEs),
ySymbols = names(MEs),
colorLabels = FALSE,
colors = blueWhiteRed(50),
textMatrix = textMatrix,
setStdMargins = FALSE,
cex.text = 0.5,
cex.lab = 0.8,
zlim = c(-1,1),
main = paste("Module-trait relationships"))
通过模块与各种表型的相关系数，可以很清楚的挑选自己感兴趣的模块进行下游分析了。直接看图，找找看有没有你感兴趣的模块~
上图中，最左侧的颜色块代表模块，最右侧的颜色条代表相关性范围。中间部分的热图中，颜色越深相关性越高，红色表示正相关，绿色表示负相关；每个单元格中的数字表示相关性和显著性。
如上图，brown模块与age_at_initial_pathologic_diagnosis性状表现为正相关且相关性最高。
此时，我们可以选择相关性最高的brown模块作为关键模块。一般，我们会按相关性的绝对值筛选最相关模块，即负相关模块也应该考虑在内。需要注意的是，grey模块中包含了所有未参与聚类的基因，因此是无效模块，不应用于后续分析。
以下引用小鬼的WGCNA图文详解(四)-模块与性状相关性 图 - 简书 (jianshu.com)
Figure 1: Module-trait associations. Each row corresponds to a module eigengene, column to a trait. Each cell contains the corresponding correlation and p-value. The table is color-coded by correlatio according to the color legend.
这张图有这几个部分：
1，横坐标：表型性状（trait）。那么，根据表型性状是连续型变量和分类变量如何数值化？
2，纵坐标：对应模块，用每个模块的eigengene来表示这个模块。那么eigengene又是什么，怎么理解这个eigengene呢？
3，图中的小格子：其中的数值代表什么？
4，每个性状与模块之间的相关性计算是否独立的：即表型放在一起分析和分开单独分析是否有不同？
1\横坐标：表型性状（trait）
教程中的数据如上图，行代表样本，列代表性状重量weight（g），长度length（cm）等。总共有134个样本，26个性状。
这里我截取了一小部分进行展示，一列代表一个性状。重量和长度都是连续性变量，直接用就好。分类变量如男女，可以男1，女0进行数值化。
2\纵坐标：对应模块的eigengene
Eigengene，即每个模块的第一主成分。
这里小编自己的理解就是这个模块有134个样本，n个基因，然后用pca主成分分析对这n个基因进行降维取其第一主成分作为这个模块的特征。
3，图中的小格子中的数值代表什么？
官方这里采用的pearson计算方法。
核心代码：moduleTraitCor = cor(MEs, datTraits, use = "p")
嗯，用的cor函数。默认的method=”pearson”。那么这里图中的小格子中的数值就代表每个性状和每个模块的特征值之间两两计算的相关性值以及对应的pvalue。
颜色表示红色越深，越正相关；绿色越深，越负相关。并且，在这里你如果理解了这个相关性的计算，那么应该就可以理解不同性状之间与模块计算相关性时是相互独立的。
模块与模块的相关性
通过聚类得到不同的模块之后，那么模块与模块之间是否也存在一定的关系呢？模块与性状之间的关系呢？
# 棕色模块与肿瘤的关系（聚类图和热图）
### 使用plotEigengeneNetworks函数再次验证棕色模块与肿瘤是否有关联。
MEs <- moduleEigengenes(datExpr, moduleColors)$eigengenes
# 这一行代码重新计算了基因表达矩阵 datExpr 的模块特征（Module Eigengenes）。
# 模块特征是每个基因模块的代表性表达模式。
# 选择感兴趣的性状：根据刚刚的相关性的图
samples <- as.data.frame(datTraits$age_at_initial_pathologic_diagnosis)
names(samples) <- "samples"
# 这一行代码将权重加到现有的模块特征中。cbind(MEs, samples) 将模块特征和你选择的性状列按列合并，
# 然后 orderMEs 函数对这些数据进行排序，生成一个新的矩阵 MET。
MET <- orderMEs(cbind(MEs, samples))
# Plot the relationships among the eigengenes and the trait
sizeGrWindow(5,9)
par(cex = 0.9)
plotEigengeneNetworks(MET, "", marDendro = c(0,4,1,2), marHeatmap = c(4,5,1,2),
cex.lab = 0.8, xLabelsAngle = 90)
### 这一段代码用于绘制模块特征与性状之间的关系。plotEigengeneNetworks 函数绘制了模块特征的网络图，其中模块特征之间的关系以及它们与选择的性状之间的关系都可以在图中看到。
图a衡量相似性用的是cor(EI;EJ)，即模块I与模块J的eigengene值的相关性，然后用1-cor(EI;EJ)来表示两个模块之间的距离（非相似度），表现在聚类树的高度（y轴的值），高度越短，两个模块越相似（负相关）。 图b使用的刻度是 (1 + cor(EI; EJ))/2，这个值越大两个模块越相似（正相关）。
这是参考资料及结论：小鬼的WGCNA图文详解(七)-eigengene network的可视化 - 简书 (jianshu.com)
模块与模块之间高度相关可以形成meta-modules，meta-modules定义的条件是这个超级模块内模块与模块之间的相关性至少要达到0.5。 从这个图可以看出，这里有三个meta-module，如red, brown and blue模块形成了一个meta-modules，它们内部之间的相关性比它们与性状weight体重的相关性还要大。
推广到自己的图里：
上图，由模块和性状age_at_initial_pathologic_diagnosis构成的聚类树，下图，对应上图的热图。两幅图表示的是一个意思，只不过是展示形式不一样这幅图表示模块与模块，模块与性状性状age_at_initial_pathologic_diagnosis之间的关系。同样的，yellow\green\red也可以形成一个meta-module
棕色模块内部GS与MM的相关性
小鬼的WGCNA图文详解(五)-GS 与MM的关联 - 简书 (jianshu.com)
TCGA数据挖掘实战（三）——WGCNA②
图剖成这几个部分：
1.横坐标：Module Membership in brown module，翻译过来就是：在棕色模块中的模块成员。模块成员是什么鬼？
2.纵坐标：Gene Significance (GS) for weight，关于体重这个性状的基因显著性。基因显著性又是什么鬼？
3.这个图是为了说明啥呢？
1.横坐标：MM
模块的eigengene和基因表达谱之间的相关性。
说的是啥呢，其实就是所有基因表达谱与这个模块的eigengene的相关性（cor）。最后是一个具有所有用来做WGCNA分析基因数长的向量，每一个值代表这个基因与模块之间的关系。如果这个值的绝对值接近0，那么这个基因就不是这个模块中的一部分，如果这个值的绝对值接近1，那么这个基因就与这个模块高度相关。
2.纵坐标：GS
GS为：基因和表型性状比如体重之间的相关性的绝对值。
总的来说，就是为了将表型特征信息与共表达网络联合起来，比如体重与哪个模块高度相关。详细一点专业一点就是：每一个基因的表达值与表型性状之间的相关性的绝对值。0表示这个基因与这个性状不相关，1表示高度相关。如果一个模块中的基因都有这个性状高度相关，那么这个模块也就与这个性状高度相关。
3、图中的每一个点
图中的每一个点代表一个基因，应该有3600个点。横坐标值表示基因与模块的相关性，纵坐标值表示基因与表型性状的相关性，这里可以看出与性状高度显著相关的基因往往是与这个性状显著相关的模块中的重要元素。
大家可以去验证一下自己的结果，如果一个性状与模块显著相关，那么这里GS与MM也会显著相关。
###棕色模块中肿瘤的GS与MM的相关性（散点图）
#我们将【基因显著性】 (Gene Significance, GS) 定义为基因与性状之间相关性的（绝对值），
#以此量化单个基因与我们感兴趣的性状（权重）之间的关联。
#对于每个模块，我们还定义了一个定量指标【模块成员】(module membership, MM），即模块特征基因与基因表达谱的相关性。
# 这样我们就可以量化阵列上所有基因与每个模块的相似性。
samples <- as.data.frame(datTraits$age_at_initial_pathologic_diagnosis)
names(samples) <- "samples"
# names (colors) of the modules
modNames <- substring(names(MEs), 3)
geneModuleMembership <- as.data.frame(cor(datExpr, MEs, use = "p"));
MMPvalue <- as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples))
names(geneModuleMembership) <- paste("MM", modNames, sep="")
names(MMPvalue) <- paste("p.MM", modNames, sep="")
geneTraitSignificance <- as.data.frame(cor(datExpr, samples, use = "p"))
GSPvalue <- as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples))
names(geneTraitSignificance) <- paste("GS.", names(samples), sep="")
names(GSPvalue) <- paste("p.GS.", names(samples), sep="")
##注意：其中前两行用来选择你感兴趣的性状，记得修改成对应的数据。
##针对感兴趣的模块，将GS与MM值可视化，代码如下：
module <- "brown"
column <- match(module, modNames);
moduleGenes <- moduleColors==module;
sizeGrWindow(7, 7);
par(mfrow = c(1,1));
verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
abs(geneTraitSignificance[moduleGenes, 1]),
xlab = paste("Module Membership in", module, "module"),
ylab = "Gene significance for sample type",
main = paste("Module membership vs. gene significance\n"),
cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
## 图中的每一个点代表一个基因。显然，图中GS 和 MM 高度相关，
## 说明与肿瘤高度显著相关的基因往往也是棕色模块中最重要（核心）的元素。
# 保存与肿瘤最相关的基因
names(datExpr)
table(moduleGenes)
sel_genes <- names(datExpr)[moduleColors=="green"]
write.table(sel_genes, file = "module_gene_name.txt", quote = F,
row.names = F,col.names = F)
write.table(sel_genes, file = "module_gene_name_11.24.txt", quote = F,
row.names = F,col.names = F)
怎么看？
关联分析：
代码中的verboseScatterplot绘制了模块成员（MM）与基因与肿瘤性状（GS）之间的散点图。每个点代表一个基因，散点图展示了模块成员与基因与肿瘤性状之间的相关性。图中的高度相关性表明与肿瘤高度显著相关的基因往往也是该模块中最重要的核心元素。
看看基因与模块的相关性（Module Membership, MM）和基因与性状的相关性（Gene Significance, GS）之间是否有某种关联。MM和GS呈正相关，说明这些与性状高度相关的基因，在关键模块中也扮演着举足轻重的角色。
相关论文的补充：Identification of Important Modules and Biomarkers in Breast Cancer Based on WGCNA （OncoTargets and Therapy 影响因子4.0）
Introduction: Breast cancer (BRCA) has the highest incidence among female malignancies, and the prognosis for these patients remains poor.
Materials and Methods: In this study, core modules and central genes related to BRCA were identified through a weighted gene co-expression network analysis (WGCNA). Gene expression profiles and clinical data of GSE25066 were obtained from the Gene Expression Omnibus (GEO) database. The result was validated with RNA-seq data from The Cancer Genome Atlas (TCGA) and Oncomine database. The top 30 key module genes with the highest intramodule connectivity were selected as the core genes (R2 = 0.40).
Results: According to TCGA and Oncomine datasets, seven genes were selected as candidate hub genes. Following further experimental verification, four hub genes (FAM171A1, NDFIP1, SKP1, and REEP5) were retained. Conclusion: We identified four hub genes as candidate biomarkers for BRCA. These hub genes may provide a theoretical basis for targeted therapy against BRCA.
技术路线：
共表达网络构建
首先，我们评估了3137个变异基因以测试其可用性，并使用R包' WGCNA '构建基因共表达网络。随后，我们构建了一个邻接矩阵来描述节点之间的关联强度。邻接矩阵的计算公式为： 在这个公式中，i和j代表两个不同的基因，Xi和Xj是它们各自的表达值。sij代表皮尔逊相关系数，aij代表两个基因之间相关性的强弱。在本研究中，我们选择了软阈值β = 10 (无标度R2 = 0.98)。随后，我们将邻接矩阵转化为拓扑重叠矩阵( TOM )。 TOM矩阵是一种通过比较两个节点与其他节点的加权相关性来定量描述节点间相似性的方法。 接下来，我们进行了层次聚类以识别模块，每个模块至少包含30个基因( min模体尺寸= 30 )。最后，计算特征基因，对模块进行层次聚类，合并相似模块( abline = 0.25)。
文章结果：数据预处理后得到GSE25066的基因表达矩阵( 12548个基因)。随后，我们选择方差大于所有四分位数的基因进行进一步分析( 3137个基因)。此外，在排除临床信息不完整的患者后，提取临床数据( 483例样本)。
我们使用皮尔逊相关系数对GSE25066中的样本进行聚类。去除离群点后，我们绘制了样本聚类树(图2 )。我们将软阈值设置为10 ( R2 = 0 . 98)来构建无标度网络。
“Next, we built the adjacency matrix and constructed the topological overlap matrix (Figure 3). Finally, nine modules were identified based on average hierarchical clustering and dynamic tree clipping (Figure 4A). The blue module was highly related to pathological grades; thus this module was selected as a clinically important module for further analysis (Figure 4B).”
接下来，我们构建了邻接矩阵进行拓扑构造重叠矩阵(图3 )。最后，基于平均层次聚类和动态树裁剪确定了9个模块(图4A )。
蓝色模块与病理分级高度相关；因此选择该模块作为临床重要模块进行进一步分析(图4B )。
（可以！！这跟我前面的意思实际上是一样的）
图4识别与乳腺癌临床特征相关的模块。
( A )所有差异表达基因的树状图基于相异性度量( 1-TOM )进行聚类。颜色波段显示了自动单块分析得到的结果。
( B )模块特征基因与BRCA临床性状相关性热图。我们选择了MEblue级别的block进行后续分析。
具有临床意义的模块识别
共表达模块是具有高度拓扑重叠相似性的基因集合。同一模块中的基因往往具有较高的共表达程度。在本研究中，我们使用了两种方法来识别与临床性状相关的重要模块。模块特征基因( ME )代表模块的第一主成分，用于描述模块在每个样本中的表达模式。模块隶属度( MM )是指基因与模块特征基因之间的相关系数，用来描述基因属于某个模块的可靠性。最后，我们计算模块与临床数据之间的相关性来识别显著的临床模块。
“Identification and Validation of Hub Genes Hub基因的识别与验证
The intramodule connectivity of a gene is equal to the sum of the degree of correlation between this gene and other genes in that module. The top 30 genes with the highest intramodule connectivity were selected as hub genes. Next, we used the Oncomine database (https://www.oncomine.org/resource/ main.html) to verify the expression, and the Kaplan-Meier plotter website (http://kmplot.com/analysis)forverification of the survival analysis to ensure the reliability of the results.”一个基因的模块内连接性等于该基因与该模块中其他基因的关联程度之和。选取模块内连通性最高的前30个基因作为hub基因。接下来，我们使用Oncomine数据库进行表达验证，并使用Kaplan - Meier绘图机网站进行生存分析，以确保结果的可靠性
参考资料汇总：
https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/Tutorials/A General Framework for Weighted Gene Co-Expression Network Analysis, Stat Appl Genet Mol Biol. 2005;4:Article17. Epub 2005 Aug 12WGCNA: an R package for weighted correlation network analysis.BMC Bioinformatics. 2008 Dec 29;9:559. doi: 10.1186/1471-2105-9-559. 小鬼的WGCNA图文详解(二)--软阈值选择 - 简书 (jianshu.com) https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/Tutorials/ A General Framework for Weighted Gene Co-Expression Network Analysis, Stat Appl Genet Mol Biol. 2005;4:Article17. Epub 2005 Aug 12
相关阅读
发表评论