数据读入、数据预处理、样本聚类、软阈值筛选、网络构建、可视化基因网络、导出网络图、关联表型数据、散点图相关性分析、分步法分析。
整体代码
library(WGCNA)
library(stringr)
setwd('E:\\AAA大创\\代码注释版\\暑假练习\\WGCNA')
WGCNA.fpkm = read.table("GSE17536.txt",header=T, comment.char = "",sep = '\t',check.names=F)
#行基因列样本，左上角（1,1）为Sample
#将文件的第一行作为列名。
#禁用注释行识别,默认情况下，R会忽略以#开头的行（如# This is a comment）。
#此参数设为空字符串后，所有行都会被当作数据读取（即使包含#）。
#禁止R自动修正列名,避免sample-1 -> sample.1#以下操作为行名作为第一列读入时进行的操作
names(WGCNA.fpkm)
datExpr0 = as.data.frame(t(WGCNA.fpkm[,-1]))
names(datExpr0) = WGCNA.fpkm$Sample;
#names改变列名
rownames(datExpr0) = names(WGCNA.fpkm[,-1])##为了训练方便减小数据规模
data1=FSbyVar(datExpr0, cut.type="topk",value=1000)
#不能用，在library(CancerSubtypes)之后var函数被重定义导致无法使用，
#连带着下面这一句也不能用，找错误找了一个多小时。。。
geneVar <- apply(datExpr0, 2, stats::var, na.rm = TRUE)
top1000 <- names(sort(geneVar, decreasing = TRUE))[1:1000]
datExpr0 <- datExpr0[, top1000]
#后续训练中所用特征数据均使用排序后前几个基因进行模拟#检查样本和基因的质量
gsg = goodSamplesGenes(datExpr0, verbose = 3)
gsg$allOK
#verbose = 3 表示输出详细的信息，包括基因和样本的质量检查。
#gsg$allOK 用来检查所有的样本和基因是否都通过了质量控制。
#如果 allOK 为 TRUE，表示数据中没有需要删除的基因或样本。#进行质量不合格数据的清除
if (!gsg$allOK) { if (sum(!gsg$goodGenes)>0) printFlush(paste("Removing genes:", paste(names(datExpr0)[!gsg$goodGenes], collapse = ", "))) if (sum(!gsg$goodSamples)>0) printFlush(paste("Removing samples:", paste(rownames(datExpr0)[!gsg$goodSamples], collapse = ", "))) datExpr0 = datExpr0[gsg$goodSamples, gsg$goodGenes]
}#计算平均表达量，过滤数据
meanFPKM = 0.5 ####--过滤标准，可以修改
#选择表达量FPKM大于0.5的基因
n=nrow(datExpr0)
datExpr0[n+1,]=apply(datExpr0[c(1:nrow(datExpr0)),],2,mean)
#2表示按列操作
datExpr0=datExpr0[1:n,datExpr0[n+1,] > meanFPKM]
#去除过小的基因#以下为保存过滤完成的数据，一般不需要运行
filtered_fpkm=t(datExpr0)
filtered_fpkm=data.frame(rownames(filtered_fpkm),filtered_fpkm)
names(filtered_fpkm)[1]="sample"
head(filtered_fpkm)
write.table(filtered_fpkm, file="mRNA.symbol.uniq.filter.txt", row.names=F, col.names=T,quote=FALSE,sep="\t")
#不对列名进行引号引用#############Sample cluster###########
sampleTree = hclust(dist(datExpr0), method = "average")
#dist(datExpr0)：计算样本之间的距离矩阵。
#hclust()使用层次聚类算法来对样本进行聚类。它的输入是样本之间的距离矩阵。
#method = "average" 选择了平均连接法，即在合并聚类时计算两个簇之间的平均距离。
pdf(file = "1.sampleClustering.pdf", width = 15, height = 8)
par(cex = 0.6)
#文字大小0.6
par(mar = c(0,6,6,0))
#图形下左上右的边距
plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 2, cex.axis = 1.5, cex.main = 2)
#main设置标题，sub设置副标题，xlab设置x轴标签，cex.lab设置两轴标签大小
#cex.axis设置坐标轴刻度标签大小 cex.main设置主标题大小
### Plot a line to show the cut
#abline(h = 180, col = "red")##剪切高度不确定，故无红线
#绘制水平红线，表示聚类的切割高度
dev.off()#添加性状数据
allTraits = read.table("type.txt",row.names=1,header=T,comment.char = "",check.names=F)
#allTraits = datExpr0[,1:3]
#这里直接采用方差最大的几个基因作为性状进行练习
dim(allTraits)
names(allTraits)
head(allTraits)
#匹配表达量和性状数据
fpkmSamples = rownames(datExpr0)
traitSamples =rownames(allTraits)
traitRows = match(fpkmSamples, traitSamples)
datTraits = allTraits[traitRows,]
rownames(datTraits)
collectGarbage() #清理内存
colorTraits <- numbers2colors(datTraits,colors = gray.colors(100, start = 0.9, end = 0.1),signed = FALSE
#FALSE处理全正数数据，TRUE处理正负数数据
)
#blueWhiteRed <- colorRampPalette(c("blue", "white", "red"))(100)
#colorTraits <- numbers2colors(datTraits, colors = blueWhiteRed, signed = FALSE)
#自定义颜色
pdf(file="2.Sample_dendrogram_and_trait_heatmap.pdf",width=20,height=12)
plotDendroAndColors(sampleTree, colors = colorTraits,groupLabels = names(datTraits), main = "Sample dendrogram and trait heatmap",cex.colorLabels = 1.5, cex.dendroLabels = 1, cex.rowText = 2)
dev.off() #软阈值筛选
powers = c(c(1:10), seq(from = 12, to=30, by=2))
sft = pickSoftThreshold(datExpr0, powerVector=powers,networkType=networkType, verbose=5)
#默认unsigned, verbose表示输出的详细程度，5最详细
par(mfrow = c(1,2))
#设置图形布局，此处为一行两列两张图
cex1 = 0.9
#文本标签大小
# 横轴是Soft threshold (power)，纵轴是无标度网络的评估参数，数值越高，
# 网络越符合无标度特征 (non-scale)
pdf(file="3.power_select.pdf",width=20,height=12)
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],#sft$fitIndices[,1]：X轴数据，表示测试的软阈值(power)值#-sign(sft$fitIndices[,3])*sft$fitIndices[,2]：Y轴数据，计算无标度拓扑拟合度#sft$fitIndices[,2]：R²值（拟合优度）#sft$fitIndices[,3]：斜率（slope）#-sign(...)：确保结果为正值（因为斜率通常为负）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],#sft$fitIndices[,5]：Y轴数据，表示平均连通性(mean connectivity)#即网络中所有基因连接度的平均值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")
dev.off()power = sft$powerEstimate
power#经验power值
# WGCNA官方建议测试power范围通常是1-30（无向）或5-30（有向）
#若出现不寻常结果，数据没有问题的话可尝试经验值
type = 'unsigned'
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)) ))
}#网络构建
##一步法网络构建：One-step network construction and module detection##
# power: 上一步计算的软阈值
# maxBlockSize: 计算机能处理的最大模块的基因数量 (默认5000)；
# 4G内存电脑可处理8000-10000个，16G内存电脑可以处理2万个，32G内存电脑可
# 以处理3万个
# 计算资源允许的情况下最好放在一个block里面。
# TOMType：拓扑重叠矩阵类型，通常为 "unsigned"（无符号相关）、
# "signed"（有符号相关）或 "signed hybrid"
# minModuleSize：最小模块大小
# 总基因数 < 5,000：较小值（15-30）
# 总基因数 5,000-20,000：中等值（30-50）
# 总基因数 > 20,000：较大值（50-100）
# reassignThreshold = 0，模块重分配阈值
# mergeCutHeight = 0.25 模块合并的切割高度
# 层级聚类中模块特征向量距离小于此值（0.25）的模块会被合并
# 基因与模块的相关性高于此阈值时会被重新分配到该模块（0 表示不重新分配）
# numericLabels: 返回数字而不是颜色作为模块的名字，后面可以再转换为颜色
# pamRespectsDendro = FALSE PAM 阶段是否尊重层级结构
# 通常设为 FALSE 以加快计算，避免 PAM 步骤覆盖初始聚类
# corType: pearson or bicor
# maxPOutliers: 设定每个基因对进行相关性计算时，允许修剪的离群样本比例上限
# 当计算两个基因的相关性时，会去除表达量分布中极端值样本
# 标准分析 0.05 默认值，平衡稳健性与信息保留
# 数据质量高/离群值少 0.01-0.03 更保守地保留样本
# 数据噪声大/离群值多 0.08-0.10 更激进地去除离群值（需谨慎，避免过度修剪）
# 禁用离群值过滤 0 不进行任何离群值修剪
# loadTOMS：加载已保存的TOM矩阵
# saveTOMs：最耗费时间的计算，存储起来，供后续使用
# verbose = 3输出中等详细信息
corType = 'pearson'
TOMType = 'unsigned'
networkType = 'unsigned'
net = blockwiseModules(datExpr0, power = power, maxBlockSize = 5000,TOMType = TOMType, minModuleSize = 30,reassignThreshold = 0, mergeCutHeight = 0.25,numericLabels = TRUE, pamRespectsDendro = FALSE,saveTOMs=TRUE, corType = corType,maxPOutliers=0.05, loadTOMs=TRUE,saveTOMFileBase = paste0('exprMat', ".tom"),verbose = 3)table(net$colors)
# 0(grey)表示未被分入任何模块#层级聚类树展示各个模块
# Convert labels to colors for plotting
moduleLabels = net$colors
moduleColors = labels2colors(moduleLabels)
# Plot the dendrogram and the module colors underneath
# 如果对结果不满意，还可以recutBlockwiseTrees，节省计算时间
pdf(file="4.cluster_tree.pdf",width=20,height=12)
plotDendroAndColors(net$dendrograms[[1]], #层次聚类结果moduleColors[net$blockGenes[[1]]],#为每个基因分配颜色标签"Module colors",#颜色条标签dendroLabels = FALSE, hang = 0.03,#是否显示基因名称，树状图叶子末端下垂高度#0：叶子末端对齐底部#0.03（推荐）：轻微下垂，增强分支可读性#0.1：明显下垂（适用于深分支结构）addGuide = TRUE, guideHang = 0.05)#是否添加参考线以及参考线下降高度
dev.off()#绘制模块之间的相关性图
# module eigengene, 可以绘制线图，作为每个模块的基因表达趋势的展示
MEs = net$MEs
#提起网络计算得到的模块特征基因，即每个模块的第一个主成分
### 不需要重新计算，改下列名字就好
### 官方教程是重新计算的，起始可以不用这么麻烦
MEs_col = MEs
colnames(MEs_col) = paste0("ME", labels2colors(as.numeric(str_replace_all(colnames(MEs),"ME",""))))
#第二行表示去掉行名中的ME，并转换为数字格式
#labels2colors将模块编号转换为颜色名MEs_col = orderMEs(MEs_col)
#"用于根据模块之间的相关性对模块顺序进行重新排列，
#使热图中的模块排序更合理、美观（类似聚类排序）"# 根据基因间表达量进行聚类所得到的各模块间的相关性图
# marDendro/marHeatmap 设置下、左、上、右的边距
#绘制模块特征基因相关性热图 + 聚类树
pdf(file="5.Eigengene_adjacency_heatmap.pdf",width=20,height=12)
plotEigengeneNetworks(MEs_col, "Eigengene adjacency heatmap", #主标题marDendro = c(3,3,2,4), #聚类树图部分的边距，下左上右marHeatmap = c(3,4,2,2), #热图部分边距plotDendrograms = T, #显示模块间聚类树xLabelsAngle = 90) #X轴文字倾斜## 如果有表型数据，也可以跟ME数据放一起，一起出图
MEs_colpheno = orderMEs(cbind(MEs_col, allTraits))
plotEigengeneNetworks(MEs_colpheno, "Eigengene adjacency heatmap",marDendro = c(3,3,2,4),marHeatmap = c(5,5,3,2), plotDendrograms = T,xLabelsAngle = 90)
dev.off()# 可视化基因网络（Topological Overlap Matrix, TOM）热图
# 如果是一步计算所有基因，或者 blocksize >= 总基因数，
# WGCNA 会把 TOM 存储在 net$TOMFiles 中
# 在这种情况下，可以直接加载已计算好的 TOM 文件，避免重复计算
load(net$TOMFiles[1], verbose = TRUE)
# 从 WGCNA 分析结果中读取第一个 TOM 文件，
# TOM 是一个基因之间的拓扑重叠矩阵（Topological Overlap Matrix）
# verbose=TRUE 表示显示加载过程中的信
# 加载后会得到一个变量名叫 TOM 的对象# 如果是分块计算且 blocksize 较小，则需要用下面这句重新计算（非常慢）
TOM = TOMsimilarityFromExpr(datExpr0, power=power, corType=corType, networkType=networkType)TOM <- as.matrix(TOM)
# 将 TOM 转换为标准的矩阵格式，方便后续处理dissTOM = 1 - TOM
# 构建 "距离矩阵"，即拓扑不相似度（dissimilarity）
# TOM 越高代表基因越相关，因此用 1 - TOM 得到“不相似度”
# 该矩阵将在热图中用于可视化基因之间的“距离”# 为了增强中等强度的连接在热图中的可视化效果，对 dissTOM 做幂次变换
# 越大的 power 会使差异更加明显。这里的 ^7 是经验选择，可根据需要调整
plotTOM = dissTOM^7
#7用于增强图像对比度，为经验值，和前面power关系不大
#可尝试不同的幂次看看热图差异，如果感觉差异不明显，适当加大；如果太夸张（太黑或太白），可以调小。diag(plotTOM) = NA
# 将矩阵对角线设为 NA，使得热图中对角线不会影响颜色分布（避免太亮或太暗）
# 因为基因与自身的距离是 0，这些点通常没有分析意义# 使用 TOMplot 函数绘制热图
# - plotTOM：用于绘图的 dissimilarity 矩阵（已做幂次变换）
# - net$dendrograms：WGCNA 中每个模块的聚类树（行列聚类使用）
# - moduleColors：每个基因对应的模块颜色，用于在热图侧边标注模块信息
# - main：图的标题
pdf(file="6.Network_heatmap_plot.pdf",width=20,height=12)
TOMplot(plotTOM, net$dendrograms, moduleColors[net$blockGenes[[1]]],main = "Network heatmap plot, all genes")
# 注意：这一步很耗时，因为需要对矩阵的行和列都做层级聚类来排序
dev.off()##导出可用于Cytoscape的网络图
probes = colnames(datExpr0)
dimnames(TOM) <- list(probes, probes)
#给TOM矩阵行列命名# 导出网络文件到点文件和边文件，Cytoscape可进行读取，
# 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)
# weighted 是否输出加权网络，设置为 TRUE 表示每条边包含连接权重，FALSE 表示只输出连接存在与否（0或1）。
# threshold 权重阈值，只输出权重大于该阈值的边。设置为 0 表示输出所有边（推荐之后用 Cytoscape 中的筛选器进一步处理）。
# nodeNames 基因名（节点名称），包含每个基因的名称或ID（如 rownames(datExpr0)）。
# nodeAttr 节点属性（模块颜色），通常为 moduleColors，表示每个基因所在模块的颜色，Cytoscape 中可用于着色或分类。##颜色分组和特征关联
trait <- "Traits.txt"
# 读入表型数据，不是必须的
# 行样本列基因
if(trait != "") {traitData <- read.table(file=trait, sep='\t', header=T, row.names=1,check.names=FALSE, comment='',quote="")sampleName = rownames(datExpr0)traitData = traitData[match(sampleName, rownames(traitData)), ]
}### 模块与表型数据关联
if (corType=="pearson") {modTraitCor = cor(MEs_col, traitData, use = "p")
# cor计算模块和表型之间的相关性矩阵，分别传入模块和表型的数据modTraitP = corPvalueStudent(modTraitCor, nSamples = dim(MEs_col)[1])
# 基于样本数计算相关性P值
} else {modTraitCorP = bicorAndPvalue(MEs_col, traitData, robustY=robustY)
# 计算两个矩阵之间每一对列的 双权重相关系数(bicor)以及对应的 P 值
# robustY意为对y中变量进行鲁棒处理，以减少极端值的影响，也有robustX参数modTraitCor = modTraitCorP$bicormodTraitP = modTraitCorP$p
}# signif表示保留几位小数
textMatrix = paste(signif(modTraitCor, 2), "\n(", signif(modTraitP, 1), ")", sep = "")
# 整理相关性值和P值格式以方便展示
dim(textMatrix) = dim(modTraitCor)
# 将字符串数组转换为二维形式
pdf(file="7.Module-trait_relationships.pdf",width=20,height=12)
labeledHeatmap(Matrix = modTraitCor, xLabels = colnames(traitData),yLabels = colnames(MEs_col),cex.lab = 0.5, #字体缩放倍数ySymbols = colnames(MEs_col), #调整图中y轴标签#默认等于yLables，后者用于告诉函数该行是哪个模块，前者用于变更显示名称colorLabels = FALSE, #不给标签染色colors = blueWhiteRed(50), #50个颜色等级textMatrix = textMatrix, #格子内部内容setStdMargins = FALSE, #不使用标准图形边距cex.text = 0.5, #格子中文本大小 zlim = c(-1,1), #设置颜色取值范围，防止红蓝分配不均main = paste("Module-trait relationships")) #图标题
dev.off()## 图中可以看到MEyellow与HEPACAM相关
## 模块内基因与表型数据关联# 性状跟模块虽然求出了相关性，可以挑选最相关的那些模块来分析，
# 但是模块本身仍然包含非常多的基因，还需进一步的寻找最重要的基因。
# 所有的模块都可以跟基因算出相关系数，所有的连续型性状也可以跟基因的表达
# 值算出相关系数。
# 如果跟性状显著相关基因也跟某个模块显著相关，那么这些基因可能就非常重要
# 。### 计算模块与基因的相关性矩阵if (corType=="pearson") {geneModuleMembership = as.data.frame(cor(datExpr0, MEs_col, use = "p"))MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples = dim(MEs_col)[1]))
} else {geneModuleMembershipA = bicorAndPvalue(datExpr0, MEs_col, robustY=robustY)geneModuleMembership = geneModuleMembershipA$bicorMMPvalue = geneModuleMembershipA$p
}# 计算性状与基因的相关性矩阵## 只有连续型性状才能进行计算，如果是离散变量，在构建样品表时就转为0-1矩阵。if (corType=="pearson") {geneTraitCor = as.data.frame(cor(datExpr0, traitData, use = "p"))geneTraitP = as.data.frame(corPvalueStudent(as.matrix(geneTraitCor), nSamples = dim(MEs_col)[1]))
} else {geneTraitCorA = bicorAndPvalue(datExpr0, traitData, robustY=robustY)geneTraitCor = as.data.frame(geneTraitCorA$bicor)geneTraitP = as.data.frame(geneTraitCorA$p)
}# 最后把两个相关性矩阵联合起来,指定感兴趣模块进行分析
module = "yellow"
pheno = "HEPACAM2"
modNames = substring(colnames(MEs_col), 3)
# 获取关注的列
module_column = match(module, modNames)
pheno_column = match(pheno,colnames(traitData))
# 获取模块内的基因
moduleGenes = moduleColors == module#sizeGrWindow(7, 7) #设置图形窗口大小
pdf(file="8.Module_membership_vs._gene_significance.pdf",width=20,height=12)
par(mfrow = c(1,1)) #设置图像布局一行一列，绘制一张图
# 与性状高度相关的基因，也是与性状相关的模型的关键基因
# WGCNA函数，该函数用于绘制带详细注释的散点图
verboseScatterplot(abs(geneModuleMembership[moduleGenes, module_column]),#模块组成基因和模块的相关性，只关注强度不关注正负，因此abs取绝对值abs(geneTraitCor[moduleGenes, pheno_column]),#模块组成基因和指定特征的相关性xlab = paste("Module Membership in", module, "module"),ylab = paste("Gene significance for", pheno),#用paste拼接xy轴标签main = paste("Module membership vs. gene significance\n"),cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)#主标题字体1.2倍 坐标轴标签1.2倍 坐标轴刻度1.2倍 散点颜色设置为模块颜色
dev.off()##教程提供了使用分步法进行模块操作的代码，即对blockwiseModules()函数进行拆分
##二者功能等价
### 计算邻接矩阵
adjacency = adjacency(datExpr0, power = power)
#计算出基因对之间的连接强度
### 把邻接矩阵转换为拓扑重叠矩阵，以降低噪音和假相关，获得距离矩阵。
TOM = TOMsimilarity(adjacency) #拓扑重叠矩阵，比简单相关性更鲁棒
dissTOM = 1-TOM### 层级聚类计算基因之间的距离树
geneTree = hclust(as.dist(dissTOM), method = "average")### 动态剪切树得到初始模块，模块合并
minModuleSize = 30
dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM,deepSplit = 2, # 控制剪枝的灵敏度，2为中等到强剪枝pamRespectsDendro = FALSE, #不使用PAM聚类校正minClusterSize = minModuleSize)
# 将模块数字标签转变为颜色标签
dynamicColors = labels2colors(dynamicMods)### 计算每个模块的特征基因（主成分分析得到特征向量）
MEList = moduleEigengenes(datExpr0, colors = dynamicColors)
MEs = MEList$eigengenes #每个模块的第一主成分
MEDiss = 1-cor(MEs)
#cor计算模块间相关性，1-cor转换为距离矩阵
METree = hclust(as.dist(MEDiss), method = "average")
#模块特征基因进行聚类##合并表达模式相似模块
MEDissThres = 0.25 #合并阈值，合并距离小于0.25的模块merge = mergeCloseModules(datExpr0, dynamicColors, cutHeight = MEDissThres, verbose = 3)
#verbose仍控制输出的详细程度，结果储存于merge中
mergedColors = merge$colors#可能因为一步法中pam修正默认开启等原因，最后得到的结果和一步法不同
一、数据读入
Show Code
library(WGCNA)
library(stringr)
setwd('E:\\AAA大创\\代码注释版\\暑假练习\\WGCNA')
WGCNA.fpkm = read.table("GSE17536.txt",header=T, comment.char = "",sep = '\t',check.names=F)
#行基因列样本，左上角（1,1）为Sample
#将文件的第一行作为列名。
#禁用注释行识别,默认情况下，R会忽略以#开头的行（如# This is a comment）。
#此参数设为空字符串后，所有行都会被当作数据读取（即使包含#）。
#禁止R自动修正列名,避免sample-1 -> sample.1#以下操作为行名作为第一列读入时进行的操作
names(WGCNA.fpkm)
datExpr0 = as.data.frame(t(WGCNA.fpkm[,-1]))
names(datExpr0) = WGCNA.fpkm$Sample;
#names改变列名
rownames(datExpr0) = names(WGCNA.fpkm[,-1])##为了训练方便减小数据规模
data1=FSbyVar(datExpr0, cut.type="topk",value=1000)
#注意在library(CancerSubtypes)之后var函数被重定义导致该函数无法使用，
#连带着下面普通的var这一句也不能用，需要加stats，找错误找了一个多小时。。。
geneVar <- apply(datExpr0, 2, stats::var, na.rm = TRUE)
top1000 <- names(sort(geneVar, decreasing = TRUE))[1:1000]
datExpr0 <- datExpr0[, top1000]
#后续训练中所用特征数据均使用排序后前几个基因进行模拟
二、数据预处理
Show Code
#检查样本和基因的质量
gsg = goodSamplesGenes(datExpr0, verbose = 3)
gsg$allOK
#verbose = 3 表示输出详细的信息，包括基因和样本的质量检查。
#gsg$allOK 用来检查所有的样本和基因是否都通过了质量控制。
#如果 allOK 为 TRUE，表示数据中没有需要删除的基因或样本。#进行质量不合格数据的清除
if (!gsg$allOK) { if (sum(!gsg$goodGenes)>0) printFlush(paste("Removing genes:", paste(names(datExpr0)[!gsg$goodGenes], collapse = ", "))) if (sum(!gsg$goodSamples)>0) printFlush(paste("Removing samples:", paste(rownames(datExpr0)[!gsg$goodSamples], collapse = ", "))) datExpr0 = datExpr0[gsg$goodSamples, gsg$goodGenes]
}#计算平均表达量，过滤数据
meanFPKM = 0.5 ####--过滤标准，可以修改
#选择表达量FPKM大于0.5的基因
n=nrow(datExpr0)
datExpr0[n+1,]=apply(datExpr0[c(1:nrow(datExpr0)),],2,mean)
#2表示按列操作
datExpr0=datExpr0[1:n,datExpr0[n+1,] > meanFPKM]
#去除过小的基因#以下为保存过滤完成的数据，一般不需要运行
filtered_fpkm=t(datExpr0)
filtered_fpkm=data.frame(rownames(filtered_fpkm),filtered_fpkm)
names(filtered_fpkm)[1]="sample"
head(filtered_fpkm)
write.table(filtered_fpkm, file="mRNA.symbol.uniq.filter.txt", row.names=F, col.names=T,quote=FALSE,sep="\t")
#不对列名进行引号引用
三、样本聚类
Show Code
#############Sample cluster###########
sampleTree = hclust(dist(datExpr0), method = "average")
#dist(datExpr0)：计算样本之间的距离矩阵。
#hclust()使用层次聚类算法来对样本进行聚类。它的输入是样本之间的距离矩阵。
#method = "average" 选择了平均连接法，即在合并聚类时计算两个簇之间的平均距离。
pdf(file = "1.sampleClustering.pdf", width = 15, height = 8)
par(cex = 0.6)
#文字大小0.6
par(mar = c(0,6,6,0))
#图形下左上右的边距
plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 2, cex.axis = 1.5, cex.main = 2)
#main设置标题，sub设置副标题，xlab设置x轴标签，cex.lab设置两轴标签大小
#cex.axis设置坐标轴刻度标签大小 cex.main设置主标题大小
### Plot a line to show the cut
#abline(h = 180, col = "red")##剪切高度不确定，故无红线
#绘制水平红线，表示聚类的切割高度
dev.off()#添加性状数据
allTraits = read.table("type.txt",row.names=1,header=T,comment.char = "",check.names=F)
#allTraits = datExpr0[,1:3]
#这里直接采用方差最大的几个基因作为性状进行练习
dim(allTraits)
names(allTraits)
head(allTraits)
#匹配表达量和性状数据
fpkmSamples = rownames(datExpr0)
traitSamples =rownames(allTraits)
traitRows = match(fpkmSamples, traitSamples)
datTraits = allTraits[traitRows,]
rownames(datTraits)
collectGarbage() #清理内存
colorTraits <- numbers2colors(datTraits,colors = gray.colors(100, start = 0.9, end = 0.1),signed = FALSE
#FALSE处理全正数数据，TRUE处理正负数数据
)
#blueWhiteRed <- colorRampPalette(c("blue", "white", "red"))(100)
#colorTraits <- numbers2colors(datTraits, colors = blueWhiteRed, signed = FALSE)
#自定义颜色
pdf(file="2.Sample_dendrogram_and_trait_heatmap.pdf",width=20,height=12)
plotDendroAndColors(sampleTree, colors = colorTraits,groupLabels = names(datTraits), main = "Sample dendrogram and trait heatmap",cex.colorLabels = 1.5, cex.dendroLabels = 1, cex.rowText = 2)
dev.off()
\(1.sampleClustering.pdf \&\& 2.Sample\_dendrogram\_and\_trait\_heatmap.pdf\)
四、软阈值筛选
Show Code
#软阈值筛选
powers = c(c(1:10), seq(from = 12, to=30, by=2))
sft = pickSoftThreshold(datExpr0, powerVector=powers,networkType=networkType, verbose=5)
#默认unsigned, verbose表示输出的详细程度，5最详细
par(mfrow = c(1,2))
#设置图形布局，此处为一行两列两张图
cex1 = 0.9
#文本标签大小
# 横轴是Soft threshold (power)，纵轴是无标度网络的评估参数，数值越高，
# 网络越符合无标度特征 (non-scale)
pdf(file="3.power_select.pdf",width=20,height=12)
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],#sft$fitIndices[,1]：X轴数据，表示测试的软阈值(power)值#-sign(sft$fitIndices[,3])*sft$fitIndices[,2]：Y轴数据，计算无标度拓扑拟合度#sft$fitIndices[,2]：R²值（拟合优度）#sft$fitIndices[,3]：斜率（slope）#-sign(...)：确保结果为正值（因为斜率通常为负）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],#sft$fitIndices[,5]：Y轴数据，表示平均连通性(mean connectivity)#即网络中所有基因连接度的平均值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")
dev.off()power = sft$powerEstimate
power
#> power
#[1] 3#经验power值
# WGCNA官方建议测试power范围通常是1-30（无向）或5-30（有向）
#若出现不寻常结果，数据没有问题的话可尝试经验值
type = 'unsigned'
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)) ))
}
\(3.power\_select.pdf\)
五、网络构建
Show Code
#网络构建
##一步法网络构建：One-step network construction and module detection##
# power: 上一步计算的软阈值
# maxBlockSize: 计算机能处理的最大模块的基因数量 (默认5000)；
# 4G内存电脑可处理8000-10000个，16G内存电脑可以处理2万个，32G内存电脑可
# 以处理3万个
# 计算资源允许的情况下最好放在一个block里面。
# TOMType：拓扑重叠矩阵类型，通常为 "unsigned"（无符号相关）、
# "signed"（有符号相关）或 "signed hybrid"
# minModuleSize：最小模块大小
# 总基因数 < 5,000：较小值（15-30）
# 总基因数 5,000-20,000：中等值（30-50）
# 总基因数 > 20,000：较大值（50-100）
# reassignThreshold = 0，模块重分配阈值
# mergeCutHeight = 0.25 模块合并的切割高度
# 层级聚类中模块特征向量距离小于此值（0.25）的模块会被合并
# 基因与模块的相关性高于此阈值时会被重新分配到该模块（0 表示不重新分配）
# numericLabels: 返回数字而不是颜色作为模块的名字，后面可以再转换为颜色
# pamRespectsDendro = FALSE PAM 阶段是否尊重层级结构
# 通常设为 FALSE 以加快计算，避免 PAM 步骤覆盖初始聚类
# corType: pearson or bicor
# maxPOutliers: 设定每个基因对进行相关性计算时，允许修剪的离群样本比例上限
# 当计算两个基因的相关性时，会去除表达量分布中极端值样本
# 标准分析 0.05 默认值，平衡稳健性与信息保留
# 数据质量高/离群值少 0.01-0.03 更保守地保留样本
# 数据噪声大/离群值多 0.08-0.10 更激进地去除离群值（需谨慎，避免过度修剪）
# 禁用离群值过滤 0 不进行任何离群值修剪
# loadTOMS：加载已保存的TOM矩阵
# saveTOMs：最耗费时间的计算，存储起来，供后续使用
# verbose = 3输出中等详细信息
corType = 'pearson'
TOMType = 'unsigned'
networkType = 'unsigned'
net = blockwiseModules(datExpr0, power = power, maxBlockSize = 5000,TOMType = TOMType, minModuleSize = 30,reassignThreshold = 0, mergeCutHeight = 0.25,numericLabels = TRUE, pamRespectsDendro = FALSE,saveTOMs=TRUE, corType = corType,maxPOutliers=0.05, loadTOMs=TRUE,saveTOMFileBase = paste0('exprMat', ".tom"),verbose = 3)table(net$colors)
# 0(grey)表示未被分入任何模块#层级聚类树展示各个模块
# Convert labels to colors for plotting
moduleLabels = net$colors
moduleColors = labels2colors(moduleLabels)
# Plot the dendrogram and the module colors underneath
# 如果对结果不满意，还可以recutBlockwiseTrees，节省计算时间
pdf(file="4.cluster_tree.pdf",width=20,height=12)
plotDendroAndColors(net$dendrograms[[1]], #层次聚类结果moduleColors[net$blockGenes[[1]]],#为每个基因分配颜色标签"Module colors",#颜色条标签dendroLabels = FALSE, hang = 0.03,#是否显示基因名称，树状图叶子末端下垂高度#0：叶子末端对齐底部#0.03（推荐）：轻微下垂，增强分支可读性#0.1：明显下垂（适用于深分支结构）addGuide = TRUE, guideHang = 0.05)#是否添加参考线以及参考线下降高度
dev.off()#绘制模块之间的相关性图
# module eigengene, 可以绘制线图，作为每个模块的基因表达趋势的展示
MEs = net$MEs
#提起网络计算得到的模块特征基因，即每个模块的第一个主成分
### 不需要重新计算，改下列名字就好
### 官方教程是重新计算的，起始可以不用这么麻烦
MEs_col = MEs
colnames(MEs_col) = paste0("ME", labels2colors(as.numeric(str_replace_all(colnames(MEs),"ME",""))))
#第二行表示去掉行名中的ME，并转换为数字格式
#labels2colors将模块编号转换为颜色名MEs_col = orderMEs(MEs_col)
#"用于根据模块之间的相关性对模块顺序进行重新排列，
#使热图中的模块排序更合理、美观（类似聚类排序）"# 根据基因间表达量进行聚类所得到的各模块间的相关性图
# marDendro/marHeatmap 设置下、左、上、右的边距
#绘制模块特征基因相关性热图 + 聚类树
pdf(file="5.Eigengene_adjacency_heatmap.pdf",width=20,height=12)
plotEigengeneNetworks(MEs_col, "Eigengene adjacency heatmap", #主标题marDendro = c(3,3,2,4), #聚类树图部分的边距，下左上右marHeatmap = c(3,4,2,2), #热图部分边距plotDendrograms = T, #显示模块间聚类树xLabelsAngle = 90) #X轴文字倾斜## 如果有表型数据，也可以跟ME数据放一起，一起出图
MEs_colpheno = orderMEs(cbind(MEs_col, allTraits))
plotEigengeneNetworks(MEs_colpheno, "Eigengene adjacency heatmap",marDendro = c(3,3,2,4),marHeatmap = c(5,5,3,2), plotDendrograms = T,xLabelsAngle = 90)
dev.off()
\(4.cluster\_tree.pdf \&\& 5.Eigengene\_adjacency\_heatmap.pdf\)
六、可视化基因网络
Show Code
# 可视化基因网络（Topological Overlap Matrix, TOM）热图
# 如果是一步计算所有基因，或者 blocksize >= 总基因数，
# WGCNA 会把 TOM 存储在 net$TOMFiles 中
# 在这种情况下，可以直接加载已计算好的 TOM 文件，避免重复计算
load(net$TOMFiles[1], verbose = TRUE)
# 从 WGCNA 分析结果中读取第一个 TOM 文件，
# TOM 是一个基因之间的拓扑重叠矩阵（Topological Overlap Matrix）
# verbose=TRUE 表示显示加载过程中的信
# 加载后会得到一个变量名叫 TOM 的对象# 如果是分块计算且 blocksize 较小，则需要用下面这句重新计算（非常慢）
TOM = TOMsimilarityFromExpr(datExpr0, power=power, corType=corType, networkType=networkType)TOM <- as.matrix(TOM)
# 将 TOM 转换为标准的矩阵格式，方便后续处理dissTOM = 1 - TOM
# 构建 "距离矩阵"，即拓扑不相似度（dissimilarity）
# TOM 越高代表基因越相关，因此用 1 - TOM 得到“不相似度”
# 该矩阵将在热图中用于可视化基因之间的“距离”# 为了增强中等强度的连接在热图中的可视化效果，对 dissTOM 做幂次变换
# 越大的 power 会使差异更加明显。这里的 ^7 是经验选择，可根据需要调整
plotTOM = dissTOM^7
#7用于增强图像对比度，为经验值，和前面power关系不大
#可尝试不同的幂次看看热图差异，如果感觉差异不明显，适当加大；如果太夸张（太黑或太白），可以调小。diag(plotTOM) = NA
# 将矩阵对角线设为 NA，使得热图中对角线不会影响颜色分布（避免太亮或太暗）
# 因为基因与自身的距离是 0，这些点通常没有分析意义# 使用 TOMplot 函数绘制热图
# - plotTOM：用于绘图的 dissimilarity 矩阵（已做幂次变换）
# - net$dendrograms：WGCNA 中每个模块的聚类树（行列聚类使用）
# - moduleColors：每个基因对应的模块颜色，用于在热图侧边标注模块信息
# - main：图的标题
pdf(file="6.Network_heatmap_plot.pdf",width=20,height=12)
TOMplot(plotTOM, net$dendrograms, moduleColors[net$blockGenes[[1]]],main = "Network heatmap plot, all genes")
# 注意：这一步很耗时，因为需要对矩阵的行和列都做层级聚类来排序
dev.off()
\(6.Network\_heatmap\_plot.pdf\)
七、导出可用于Cytoscope的网络图
Show Code
##导出可用于Cytoscape的网络图
probes = colnames(datExpr0)
dimnames(TOM) <- list(probes, probes)
#给TOM矩阵行列命名# 导出网络文件到点文件和边文件，Cytoscape可进行读取，
# 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)
# weighted 是否输出加权网络，设置为 TRUE 表示每条边包含连接权重，FALSE 表示只输出连接存在与否（0或1）。
# threshold 权重阈值，只输出权重大于该阈值的边。设置为 0 表示输出所有边（推荐之后用 Cytoscape 中的筛选器进一步处理）。
# nodeNames 基因名（节点名称），包含每个基因的名称或ID（如 rownames(datExpr0)）。
# nodeAttr 节点属性（模块颜色），通常为 moduleColors，表示每个基因所在模块的颜色，Cytoscape 中可用于着色或分类。
导出的网络文件：
八、关联表型数据
Show Code
##颜色分组和特征关联
trait <- "Traits.txt"
# 读入表型数据，不是必须的
# 行样本列基因
if(trait != "") {traitData <- read.table(file=trait, sep='\t', header=T, row.names=1,check.names=FALSE, comment='',quote="")sampleName = rownames(datExpr0)traitData = traitData[match(sampleName, rownames(traitData)), ]
}### 模块与表型数据关联
if (corType=="pearson") {modTraitCor = cor(MEs_col, traitData, use = "p")
# cor计算模块和表型之间的相关性矩阵，分别传入模块和表型的数据modTraitP = corPvalueStudent(modTraitCor, nSamples = dim(MEs_col)[1])
# 基于样本数计算相关性P值
} else {modTraitCorP = bicorAndPvalue(MEs_col, traitData, robustY=robustY)
# 计算两个矩阵之间每一对列的 双权重相关系数(bicor)以及对应的 P 值
# robustY意为对y中变量进行鲁棒处理，以减少极端值的影响，也有robustX参数modTraitCor = modTraitCorP$bicormodTraitP = modTraitCorP$p
}# signif表示保留几位小数
textMatrix = paste(signif(modTraitCor, 2), "\n(", signif(modTraitP, 1), ")", sep = "")
# 整理相关性值和P值格式以方便展示
dim(textMatrix) = dim(modTraitCor)
# 将字符串数组转换为二维形式
pdf(file="7.Module-trait_relationships.pdf",width=20,height=12)
labeledHeatmap(Matrix = modTraitCor, xLabels = colnames(traitData),yLabels = colnames(MEs_col),cex.lab = 0.5, #字体缩放倍数ySymbols = colnames(MEs_col), #调整图中y轴标签#默认等于yLables，后者用于告诉函数该行是哪个模块，前者用于变更显示名称colorLabels = FALSE, #不给标签染色colors = blueWhiteRed(50), #50个颜色等级textMatrix = textMatrix, #格子内部内容setStdMargins = FALSE, #不使用标准图形边距cex.text = 0.5, #格子中文本大小 zlim = c(-1,1), #设置颜色取值范围，防止红蓝分配不均main = paste("Module-trait relationships")) #图标题
dev.off()
\(7.Module-trait\_relationships.pdf\)
九、散点图相关性分析
Show Code
## 图中可以看到MEyellow与HEPACAM相关
## 模块内基因与表型数据关联# 性状跟模块虽然求出了相关性，可以挑选最相关的那些模块来分析，
# 但是模块本身仍然包含非常多的基因，还需进一步的寻找最重要的基因。
# 所有的模块都可以跟基因算出相关系数，所有的连续型性状也可以跟基因的表达
# 值算出相关系数。
# 如果跟性状显著相关基因也跟某个模块显著相关，那么这些基因可能就非常重要
# 。### 计算模块与基因的相关性矩阵if (corType=="pearson") {geneModuleMembership = as.data.frame(cor(datExpr0, MEs_col, use = "p"))MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples = dim(MEs_col)[1]))
} else {geneModuleMembershipA = bicorAndPvalue(datExpr0, MEs_col, robustY=robustY)geneModuleMembership = geneModuleMembershipA$bicorMMPvalue = geneModuleMembershipA$p
}# 计算性状与基因的相关性矩阵## 只有连续型性状才能进行计算，如果是离散变量，在构建样品表时就转为0-1矩阵。if (corType=="pearson") {geneTraitCor = as.data.frame(cor(datExpr0, traitData, use = "p"))geneTraitP = as.data.frame(corPvalueStudent(as.matrix(geneTraitCor), nSamples = dim(MEs_col)[1]))
} else {geneTraitCorA = bicorAndPvalue(datExpr0, traitData, robustY=robustY)geneTraitCor = as.data.frame(geneTraitCorA$bicor)geneTraitP = as.data.frame(geneTraitCorA$p)
}# 最后把两个相关性矩阵联合起来,指定感兴趣模块进行分析
module = "yellow"
pheno = "HEPACAM2"
modNames = substring(colnames(MEs_col), 3)
# 获取关注的列
module_column = match(module, modNames)
pheno_column = match(pheno,colnames(traitData))
# 获取模块内的基因
moduleGenes = moduleColors == module#sizeGrWindow(7, 7) #设置图形窗口大小
pdf(file="8.Module_membership_vs._gene_significance.pdf",width=20,height=12)
par(mfrow = c(1,1)) #设置图像布局一行一列，绘制一张图
# 与性状高度相关的基因，也是与性状相关的模型的关键基因
# WGCNA函数，该函数用于绘制带详细注释的散点图
verboseScatterplot(abs(geneModuleMembership[moduleGenes, module_column]),#模块组成基因和模块的相关性，只关注强度不关注正负，因此abs取绝对值abs(geneTraitCor[moduleGenes, pheno_column]),#模块组成基因和指定特征的相关性xlab = paste("Module Membership in", module, "module"),ylab = paste("Gene significance for", pheno),#用paste拼接xy轴标签main = paste("Module membership vs. gene significance\n"),cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)#主标题字体1.2倍 坐标轴标签1.2倍 坐标轴刻度1.2倍 散点颜色设置为模块颜色
dev.off()
\(8.Module\_membership_vs.\_gene\_significance.pdf\)（恰好是黄色的点，看着不大明显）
十、分步法分析
Show Code
##网络教程提供了使用分步法进行模块操作的代码，即对blockwiseModules()函数进行拆分
##二者功能等价
### 计算邻接矩阵
adjacency = adjacency(datExpr0, power = power)
#计算出基因对之间的连接强度
### 把邻接矩阵转换为拓扑重叠矩阵，以降低噪音和假相关，获得距离矩阵。
TOM = TOMsimilarity(adjacency) #拓扑重叠矩阵，比简单相关性更鲁棒
dissTOM = 1-TOM### 层级聚类计算基因之间的距离树
geneTree = hclust(as.dist(dissTOM), method = "average")### 动态剪切树得到初始模块，模块合并
minModuleSize = 30
dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM,deepSplit = 2, # 控制剪枝的灵敏度，2为中等到强剪枝pamRespectsDendro = FALSE, #不使用PAM聚类校正minClusterSize = minModuleSize)
# 将模块数字标签转变为颜色标签
dynamicColors = labels2colors(dynamicMods)### 计算每个模块的特征基因（主成分分析得到特征向量）
MEList = moduleEigengenes(datExpr0, colors = dynamicColors)
MEs = MEList$eigengenes #每个模块的第一主成分
MEDiss = 1-cor(MEs)
#cor计算模块间相关性，1-cor转换为距离矩阵
METree = hclust(as.dist(MEDiss), method = "average")
#模块特征基因进行聚类##合并表达模式相似模块
MEDissThres = 0.25 #合并阈值，合并距离小于0.25的模块merge = mergeCloseModules(datExpr0, dynamicColors, cutHeight = MEDissThres, verbose = 3)
#verbose仍控制输出的详细程度，结果储存于merge中
mergedColors = merge$colors#可能因为一步法中pam修正默认开启等原因，最后得到的结果和一步法不同