本節(jié)概覽
1. WGCNA基本概念：
定義、關(guān)鍵術(shù)語、基本流程、一些注意事項
2. WGCNA運行：
?輸入數(shù)據(jù)準備
①判斷數(shù)據(jù)質(zhì)量，繪制樣品的系統(tǒng)聚類樹
②挑選最佳閾值power
③ 構(gòu)建加權(quán)共表達網(wǎng)絡(luò)（ 一步法和分步法），識別基因模塊
④ 關(guān)聯(lián)基因模塊與表型：模塊與表型相關(guān)性熱圖、模塊與表型相關(guān)性boxplot圖、基因與模塊、表型相關(guān)性散點圖
⑤ WGCNA的標(biāo)配熱圖 ，模塊相關(guān)性展示
⑥ 對感興趣模塊的基因進行批量GO分析
⑦ 感興趣模塊繪制熱圖
⑧ 提取感興趣模塊的基因名, 導(dǎo)出基因至 VisANT 或 cytoscape作圖
WGCNA的本質(zhì)是“從基因的協(xié)同表達模式出發(fā)，把基因按‘共表達團伙’進行分組，再找出哪些團伙的活動規(guī)律與樣本的生理、病理或處理條件等表型指標(biāo)存在強關(guān)聯(lián)”。
也就是說，它不關(guān)心單個基因的表達差異，而是先通過表達相似性把基因聚成若干個“模塊”（每個模塊相當(dāng)于一個基因小團體），然后計算每個模塊的“綜合表達特征”與各種表型（比如疾病 VS 對照、連續(xù)變量用藥劑量、時間點等）的相關(guān)性。最終，那些與特定表型顯著相關(guān)的模塊，就被認為是該表型潛在的關(guān)鍵功能單元。
1. WGCNA基本概念
1.1 定義
WGCNA 就是把海量基因按照它們在整個樣本集中的表達“步調(diào)一致性”來分組，形成一個個功能“陣營”（模塊）；每個陣營可以提煉出一個“代表得分”（eigengene）或挑出一個“樞紐基因”（hub gene）來代言該陣營的行為，然后去看每個陣營的代表得分與樣本的各種臨床、生理或?qū)嶒炛笜?biāo)是否存在統(tǒng)計學(xué)上的關(guān)聯(lián)。
1.2 其他關(guān)鍵術(shù)語
關(guān)鍵術(shù)語（Connectivity、Adjacency matrix、TOM等）見：
Simulated-00-Background.pdf (ucla.edu)
一文看懂WGCNA 分析
WGCNA分析，簡單全面的最新教程
1.3 基本流程
構(gòu)建基因共表達網(wǎng)絡(luò) >> 識別基因模塊 >> 關(guān)聯(lián)基因模塊與表型 >> 研究基因模塊間關(guān)系 >> 從感興趣的基因模塊中尋找關(guān)鍵驅(qū)動基因
1.4 一些注意事項
WGCNA package: Frequently Asked Questions (ucla.edu)
- 樣本數(shù) >= 15
-
基因過濾方法
先看基因的表達絕對值（如均值、中位數(shù)）或樣本間的波動幅度（如方差、MAD），把那些本身表達量很低、或者在所有樣本里都“平平無奇”沒什么變化的基因篩掉；但不要事先拿“不同組間有沒有差異”作為過濾標(biāo)準，那是差異分析的任務(wù)，不適合用來做預(yù)處理過濾。 -
輸入數(shù)據(jù)形式
分析前，若存在批次效應(yīng)要先校正；
針對RNA-seq數(shù)據(jù)，推薦使用DESeq2包中的方差穩(wěn)定變換（varianceStabilizingTransformation），或者對已標(biāo)準化的表達量（比如FPKM、CPM）做 log2(x+1)轉(zhuǎn)化 -
經(jīng)驗軟閾值power
如果在無向網(wǎng)絡(luò)中嘗試的 power 值小于 15、有向網(wǎng)絡(luò)中 power 值小于 30 的范圍內(nèi)，所有取值都無法同時滿足兩個條件——即無標(biāo)度拓撲擬合指數(shù) R^2≥0.8且平均連接度低于 100——那么就只能依據(jù)經(jīng)驗人為選取一個 power 值。
2. WGCNA運行
主要參考修改自以下代碼：
GEO/task4-NPC at master · jmzeng1314/GEO · GitHub
手把手10分文章WGCNA復(fù)現(xiàn)：小膠質(zhì)細胞亞群在腦發(fā)育時髓鞘形成的作用
以下實戰(zhàn)數(shù)據(jù)來自于https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE154290，Supplementary file中的GSE154290_RNA.FPKM.txt.gz文件
? 輸入數(shù)據(jù)準備
讀取fpkm表達矩陣，進行l(wèi)og2(x+1)轉(zhuǎn)化；最后對矩陣進行轉(zhuǎn)置，變成行為樣品、列為基因的形式
# 清除當(dāng)前環(huán)境中的所有變量
rm(list = ls())
# 防止字符串自動轉(zhuǎn)換為因子（保持字符類型）
options(stringsAsFactors = F)
# 設(shè)置系統(tǒng)語言為英文（避免中文報錯信息亂碼）
Sys.setenv(LANGUAGE = "en")
# 加載所需R包
library(WGCNA) # 加權(quán)基因共表達網(wǎng)絡(luò)分析
library(FactoMineR) # 多元探索性數(shù)據(jù)分析（如PCA）
library(factoextra) # 可視化FactoMineR結(jié)果
library(tidyverse) # 包含ggplot2、dplyr、tidyr等數(shù)據(jù)處理和可視化工具
library(data.table) # 多核讀取大文件（此處未直接使用但預(yù)留）
# 創(chuàng)建并切換到結(jié)果輸出目錄
dir.create('8.WGCNA')
setwd('8.WGCNA')
### 啟用WGCNA多核計算（加速分析）
# 使用可用物理核心數(shù)的75%進行并行計算
enableWGCNAThreads(nThreads = 0.75 * parallel::detectCores())
################################# 0.輸入數(shù)據(jù)準備 ################################
### 讀取表達矩陣（FPKM值）
if (T) {
# 讀取壓縮的FPKM文件（數(shù)據(jù)框格式，不轉(zhuǎn)為data.table）
fpkm00 <- fread("../GSE154290_RNA.FPKM.txt.gz", data.table = F)
# 統(tǒng)計基因名重復(fù)個數(shù)（同個基因可能對應(yīng)多行）
table(duplicated(fpkm00$gene))
# 保存基因名列，并移除它（后續(xù)用于聚合）
gene <- fpkm00$gene
fpkm00 <- fpkm00[, -1]
# 按基因名合并表達量（相同基因求和，適用于多個轉(zhuǎn)錄本合并為基因水平）
fpkm0 <- aggregate(fpkm00, by = list(gene), FUN = sum)
# 將第一列（基因名）設(shè)為行名
fpkm <- column_to_rownames(fpkm0, "Group.1")
}
# 對FPKM值進行l(wèi)og2轉(zhuǎn)換（加1避免log2(0)）
data <- log2(fpkm + 1)
### 篩選MAD（絕對中位差）前5000的基因
# MAD是比方差更穩(wěn)健的離散度度量，常用于篩選變化最顯著的基因
keep_data <- data[order(apply(data, 1, mad), decreasing = T)[1:5000], ]
### 創(chuàng)建datTraits，包含樣本的分組/表型信息
datTraits <- data.frame(row.names = colnames(data), # 行名為樣本名
group = colnames(data)) # 暫將樣本名作為分組信息（需修改）
fix(datTraits) # 彈出交互式表格手動修改分組信息（這一步在腳本中需人工干預(yù)）
### 為分組添加數(shù)值編號（WGCNA要求表型為數(shù)值型）
grouptype <- data.frame(group = sort(unique(datTraits$group)), # 所有分組名稱
groupNo = 1:length(unique(datTraits$group))) # 賦予編號
# 初始化groupNo列
datTraits$groupNo <- "NA"
# 根據(jù)分組名稱匹配編號
for(i in 1:nrow(grouptype)){
datTraits[which(datTraits$group == grouptype$group[i]), 'groupNo'] <- grouptype$groupNo[i]
}
datTraits # 查看最終的表型數(shù)據(jù)框
### 轉(zhuǎn)置表達矩陣：行為樣本，列為基因（WGCNA的標(biāo)準輸入格式）
datExpr0 <- as.data.frame(t(keep_data))
此處根據(jù)細胞的多能性狀態(tài)，將它們大致劃分為3個組，并按照發(fā)育順序依次標(biāo)記為1、2、3，然后將這一分組信息添加到樣本的表型數(shù)據(jù)中。
① 判斷數(shù)據(jù)質(zhì)量，繪制樣品的系統(tǒng)聚類樹
先評估樣本與基因的數(shù)據(jù)質(zhì)量，剔除低質(zhì)量數(shù)據(jù)；再對樣本進行層次聚類并繪制聚類樹，如果某個樣本明顯偏離整體（表現(xiàn)為一個顯著的離群點），就需要將其移除；此外，還可以通過PCA圖來觀察樣本的整體分布情況。
############################## 1.判斷數(shù)據(jù)質(zhì)量 ################################
### 判斷數(shù)據(jù)質(zhì)量--缺失值
# 檢查表達矩陣中的基因和樣本是否存在過多缺失值（或無效表達值）
# goodSamplesGenes 會篩選出缺失值比例低于閾值的樣本和基因
gsg <- goodSamplesGenes(datExpr0, verbose = 3)
gsg$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 = ", ")));
# 從數(shù)據(jù)中剔除不合格的基因和樣本
datExpr0 = datExpr0[gsg$goodSamples, gsg$goodGenes]
}
# 再次檢查，確保所有數(shù)據(jù)均合格
gsg <- goodSamplesGenes(datExpr0, verbose = 3)
gsg$allOK # 此時應(yīng)為 TRUE
### 繪制樣品的系統(tǒng)聚類樹
if(T){
# 基于樣本表達譜計算歐氏距離，并用平均鏈接法進行層次聚類
sampleTree <- hclust(dist(datExpr0), method = "average")
# 設(shè)置圖形參數(shù)：下邊距0，左邊距5，上邊距2，右邊距0
par(mar = c(0,5,2,0))
# 繪制樣本聚類樹（不展示在文件中，僅屏幕預(yù)覽）
plot(sampleTree, main = "Sample clustering", sub="", xlab="", cex.lab = 2,
cex.axis = 1, cex.main = 1, cex.lab=1)
## 為每個樣本分配顏色：根據(jù)分組信息（datTraits$group）轉(zhuǎn)換為顏色條
# 將分組因子轉(zhuǎn)為數(shù)值，再映射到彩虹色上（顏色數(shù)量等于分組數(shù)）
sample_colors <- numbers2colors(as.numeric(factor(datTraits$group)),
colors = rainbow(length(table(datTraits$group))),
signed = FALSE)
## 將聚類樹和顏色條繪制在同一張圖中，便于觀察樣本聚類與表型分組的一致性
par(mar = c(1,4,3,1), cex=0.8)
pdf("step1_Sample dendrogram and trait.pdf", width = 8, height = 6)
plotDendroAndColors(sampleTree, sample_colors,
groupLabels = "trait", # 顏色條對應(yīng)的標(biāo)簽
cex.dendroLabels = 0.8,
marAll = c(1, 4, 3, 1),
cex.rowText = 0.01,
main = "Sample dendrogram and trait" )
## 如需剔除離群樣本，可以在聚類樹高度上畫一條切割線（例子中注釋掉了）
# abline(h = 23500, col = "red") # 切割高度需根據(jù)實際情況調(diào)整
dev.off() # 關(guān)閉pdf設(shè)備
}
## 若存在顯著的離群樣本，可在此步驟剔除（當(dāng)前設(shè)置為 FALSE 表示不執(zhí)行）
if(F){
# cutreeStatic 基于樹高將樣本分成若干簇，只保留第一簇（clust==1）
clust <- cutreeStatic(sampleTree, cutHeight = 23500, minSize = 10)
table(clust)
keepSamples <- (clust==1)
datExpr0 <- datExpr0[keepSamples, ]
datTraits <- datTraits[keepSamples, ]
dim(datExpr0)
}
### 判斷數(shù)據(jù)質(zhì)量 : PCA進行分組查看
# 提取分組信息（因子或字符向量）
group_list <- datTraits$group
# 對表達矩陣（行為樣本，列為基因）執(zhí)行主成分分析
# PCA 可評估樣本的總體分布，判斷批次效應(yīng)或離群樣本
dat.pca <- PCA(datExpr0, graph = F)
# 使用 factoextra 包可視化 PCA 結(jié)果：樣本點圖，按分組著色
pca <- fviz_pca_ind(dat.pca,
title = "Principal Component Analysis",
legend.title = "Groups",
geom.ind = c("point","text"), # 同時顯示點和樣本標(biāo)簽
pointsize = 2,
labelsize = 4,
repel = TRUE, # 避免標(biāo)簽重疊
col.ind = group_list, # 按分組顏色區(qū)分
axes.linetype=NA, # 不顯示坐標(biāo)軸線
mean.point=FALSE # 不顯示各組的中心點
) +
theme(legend.position = "none") + # 移除圖例（如果不需要）
coord_fixed(ratio = 1) # 固定坐標(biāo)軸比例（使形狀不變形）
pca # 屏幕顯示圖形
# 保存 PCA 圖為 PDF 文件
ggsave(pca, filename = "step1_Sample PCA analysis.pdf", width = 8, height = 8)
## 保存當(dāng)前處理后的數(shù)據(jù)，供后續(xù)分析步驟（如選擇軟閾值、構(gòu)建網(wǎng)絡(luò)）使用
datExpr <- datExpr0
nGenes <- ncol(datExpr)
nSamples <- nrow(datExpr)
save(nGenes, nSamples, datExpr, datTraits, file = "step1_input.Rdata")
下圖可看出此次聚類效果還是比較好的，沒有離群值，因此不用進行剔除樣本
② 挑選最佳閾值power
挑選標(biāo)準：R^2 > 0.8 , slope ≈ -1，選擇R^2 與 power關(guān)系圖的拐點處的power值
https://horvath.genetics.ucla.edu/html/GeneralFramework/YeastTutorialHorvath.pdf
sft$powerEstimate可以給出所設(shè)定R^2 cut-off（默認為0.85）處的power值
############################### 2.挑選最佳閾值power ###################################
# 清空環(huán)境，避免之前變量的干擾
rm(list = ls())
# 加載第一步保存的表達矩陣和表型數(shù)據(jù)（datExpr, datTraits）
load("step1_input.Rdata")
# 設(shè)置無標(biāo)度拓撲擬合指數(shù)的閾值（R^2），通常使用0.8或0.85
R.sq_cutoff = 0.8
# 執(zhí)行軟閾值冪次篩選（若條件始終為真，可省略 if(T) 直接運行）
if(T){
# 定義候選冪次（power）范圍：1~20連續(xù)，22~30以2為步長
powers <- c(seq(1,20,by = 1), seq(22,30,by = 2))
# 調(diào)用WGCNA函數(shù)評估不同冪次下的網(wǎng)絡(luò)拓撲特性
# networkType = "unsigned" 表示無符號網(wǎng)絡(luò)（不考慮正負相關(guān)性）
# RsquaredCut 為提前設(shè)定的R^2閾值，函數(shù)會嘗試尋找滿足條件的最小冪次
sft <- pickSoftThreshold(datExpr,
networkType = "unsigned",
powerVector = powers,
RsquaredCut = R.sq_cutoff,
verbose = 5) # verbose=5輸出詳細過程
# 理想的無標(biāo)度網(wǎng)絡(luò)要求R^2 > 0.8 且 平均連接度（mean connectivity）下降的斜率接近-1
# 將結(jié)果保存為PDF文件（寬16英寸，高12英寸）
pdf("step2_power-value.pdf", width = 16, height = 12)
# 設(shè)置圖形布局為一行兩列（并排兩個圖）
par(mfrow = c(1,2));
cex1 = 0.9; # 后續(xù)文本字符大小
# 左圖：無標(biāo)度拓撲擬合指數(shù)（signed R^2）隨冪次的變化
# sft$fitIndices[,1] -> 冪次；sft$fitIndices[,2] -> 無標(biāo)度R^2；sft$fitIndices[,3] -> 斜率符號
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")
# 在圖上標(biāo)注每個冪次的數(shù)值（紅色）
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
labels=powers, cex=cex1, col="red")
# 添加R^2閾值線（紅色水平線）
abline(h=R.sq_cutoff, col="red")
# 右圖：平均連接度（Mean Connectivity）隨冪次的變化
# sft$fitIndices[,5] 即為平均連接度
plot(sft$fitIndices[,1], sft$fitIndices[,5],
xlab="Soft Threshold (power)", ylab="Mean Connectivity", type="n")
text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1, col="red")
# 平均連接度參考線（通常希望連接度不要過高，例如低于100）
abline(h=100, col="red")
dev.off() # 關(guān)閉PDF設(shè)備
}
# 輸出函數(shù)自動估計的最佳冪次（若滿足R^2閾值則返回對應(yīng)冪次，否則返回NA）
sft$powerEstimate
# 此處手動設(shè)定冪次為15（可根據(jù)圖形判斷或sft$powerEstimate結(jié)果修改）
# power = sft$powerEstimate
power = 15
# 當(dāng)自動估計的冪次為NA（即沒有任何冪次滿足R^2>0.8且平均連接度<100）時，使用經(jīng)驗公式計算冪次
# 這可能由于樣品間差異過大、批次效應(yīng)或存在離群樣本導(dǎo)致無標(biāo)度擬合不理想
if(is.na(power)){
# 網(wǎng)絡(luò)類型（與pickSoftThreshold參數(shù)一致，此處為"unsigned"）
type = "unsigned"
nSamples = nrow(datExpr) # 樣本數(shù)量
# 根據(jù)樣本數(shù)量和經(jīng)驗查表給出冪次建議（unsigned網(wǎng)絡(luò)下樣本少時用較小冪次）
# 若使用signed或signed hybrid網(wǎng)絡(luò)，冪次通常需增大（參考官方文檔）
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))
)
)
}
# 保存冪次篩選結(jié)果（sft對象）和最終選定的power值，供后續(xù)網(wǎng)絡(luò)構(gòu)建使用
save(sft, power, file = "step2_power_value.Rdata")
圖中可見，拐點約在 power 取 15 處，且各參數(shù)基本達到標(biāo)準，因此選擇 power = 15。
③ 構(gòu)建加權(quán)共表達網(wǎng)絡(luò)，識別基因模塊
構(gòu)建加權(quán)共表達網(wǎng)絡(luò)時，既可以采用一步法，也可以采用分步法。通常為了簡便會優(yōu)先選擇一步法；但如果希望自由調(diào)節(jié)最終獲得的模塊數(shù)量，則分步法更加靈活。
對于一步法構(gòu)建網(wǎng)絡(luò)并識別模塊而言，核心可調(diào)參數(shù)是 minModuleSize 和 mergeCutHeight，這兩個參數(shù)取值越大，得到的模塊數(shù)量就越少。
##################### 3.一步法構(gòu)建加權(quán)共表達網(wǎng)絡(luò)，識別基因模塊 ####################
# 加載第一步保存的表達矩陣和表型數(shù)據(jù)（datExpr, datTraits）
load(file = "step1_input.Rdata")
# 加載第二步確定的軟閾值冪次（power）
load(file = "step2_power_value.Rdata")
# 執(zhí)行一步法網(wǎng)絡(luò)構(gòu)建和模塊檢測（條件始終為真，可省略 if(T)）
if(T){
net <- blockwiseModules(
datExpr, # 表達矩陣：行為樣本，列為基因
power = power, # 軟閾值冪次（由上一步確定）
# 分塊設(shè)置：將所有基因放在同一個塊中計算（若基因數(shù)過多可設(shè)置較小值如5000）
maxBlockSize = ncol(datExpr), # 此處直接設(shè)置為總基因數(shù)，表示不拆分
corType = "pearson", # 計算相關(guān)性的方法："pearson"（默認）或 "bicor"
# "bicor" 對離群點更穩(wěn)健，但計算量更大
networkType = "unsigned", # 鄰接矩陣類型："unsigned"（只考慮正相關(guān)，忽略負相關(guān)）
# 可選 "signed"（保留正負符號）、"signed hybrid"
TOMType = "unsigned", # TOM（拓撲重疊矩陣）的計算類型：
# "unsigned" 忽略正負號，通常與 unsigned 網(wǎng)絡(luò)匹配
minModuleSize = 30, # 模塊最小基因數(shù)：小于該值的模塊將被合并到其他模塊或灰色組
# 值越大，最終得到的模塊數(shù)量越少
mergeCutHeight = 0.25, # 模塊合并閾值：動態(tài)樹切割后，將相似度高于該閾值的模塊合并
# 值越大，合并越多，模塊數(shù)越少（范圍0~1，典型值0.25）
numericLabels = TRUE, # 模塊標(biāo)簽返回數(shù)字（TRUE）還是顏色（FALSE）
# 數(shù)字便于后續(xù)計算，可再用 labels2colors() 轉(zhuǎn)為顏色
saveTOMs = FALSE, # 是否保存TOM矩陣到文件（占用大量磁盤空間，常設(shè)為FALSE）
verbose = 3 # 輸出詳細日志級別：0靜默，3較詳細，5最詳細
)
# 查看每個模塊包含的基因數(shù)量（0 表示灰色模塊，即未分配任何模塊的基因）
table(net$colors)
# 參數(shù)提示（原代碼中的英文注釋，已轉(zhuǎn)化為中文補充說明）：
# power: 上一步計算的軟閾值
# maxBlockSize: 每個塊能處理的最大基因數(shù)（取決于內(nèi)存大小，16G內(nèi)存可處理約2萬個基因）
# 資源允許時最好放在一個塊內(nèi)（即 maxBlockSize >= 總基因數(shù)）
# corType：相關(guān)性計算方法；pearson（默認）或 bicor（對離群點更穩(wěn)?。? # networkType：鄰接矩陣是否考慮正負相關(guān)；unsigned（只正相關(guān)）、signed（保留符號）、signed hybrid
# TOMType：TOM矩陣是否考慮正負相關(guān)；一般與 networkType 保持一致，但因鄰接矩陣非負，所以影響不大
# numericLabels: 模塊編號為數(shù)字（而非顏色），便于程序處理，后續(xù)可用 labels2colors() 轉(zhuǎn)換
# saveTOMs：TOM計算最耗時，若需重復(fù)使用可保存為 .Rdata 文件
# mergeCutHeight: 合并相似模塊的閾值，越大模塊越少，一般 0.25 左右
# minModuleSize: 每個模塊的最小基因數(shù)，設(shè)定越大模塊越少
# 輸出模塊按基因數(shù)量降序排列，編號 1,2,3,...；灰色模塊（編號0）表示未分入任何模塊的基因
}
## 模塊可視化：使用層次聚類樹展示基因模塊劃分
if(T){
# 將數(shù)字標(biāo)簽轉(zhuǎn)換為顏色標(biāo)簽（更直觀）
moduleColors <- labels2colors(net$colors)
table(moduleColors) # 查看每種顏色模塊包含的基因數(shù)
# 繪制基因聚類樹并在下方顯示模塊顏色條
pdf("step3_genes-modules_ClusterDendrogram.pdf", width = 16, height = 12)
plotDendroAndColors(
net$dendrograms[[1]], # 基因聚類樹（第一個塊）
moduleColors[net$blockGenes[[1]]], # 對應(yīng)模塊顏色（按聚類樹葉子順序排列）
"Module colors", # 顏色條的圖例標(biāo)簽
dendroLabels = FALSE, # 不顯示樹葉標(biāo)簽（基因名太多）
hang = 0.03, # 樹枝懸掛高度
addGuide = TRUE, # 添加輔助線（分隔顏色條區(qū)域）
guideHang = 0.05 # 輔助線懸空高度
)
dev.off() # 關(guān)閉PDF設(shè)備
}
# 保存網(wǎng)絡(luò)對象和模塊顏色，供后續(xù)分析（如模塊與表型關(guān)聯(lián)、基因富集等）
save(net, moduleColors, file = "step3_genes_modules.Rdata")
- 分步法構(gòu)建加權(quán)共表達網(wǎng)絡(luò)，識別基因模塊
與一步法類似，主要調(diào)整參數(shù)為minModuleSize和MEDissThres, 數(shù)值越大模塊越少
##################### 分布法完成網(wǎng)絡(luò)構(gòu)建，一般不用 #################################
if(F){
## 構(gòu)建加權(quán)共表達網(wǎng)絡(luò)分為兩步（分布法）：
## 步驟1：計算鄰接矩陣（基因間相關(guān)性的軟閾值冪次轉(zhuǎn)換）
## 步驟2：計算拓撲重疊矩陣（TOM），并用TOM不相似度進行聚類和模塊識別
### (1) 網(wǎng)絡(luò)構(gòu)建：計算相關(guān)性 → 鄰接矩陣
# adjacency() 首先計算基因間 Pearson 相關(guān)系數(shù)，再通過 power 軟閾值將其轉(zhuǎn)換為鄰接值
# 公式：adjacency = |cor(gene_i, gene_j)|^power
# 輸出為對稱矩陣，值在[0,1]之間，表示基因?qū)Φ倪B接強度
adjacency = adjacency(datExpr, power = power)
### (2) 將鄰接矩陣轉(zhuǎn)換為拓撲重疊矩陣（TOM）
# TOM 衡量兩個基因在共表達網(wǎng)絡(luò)中的“共享鄰居”程度，能更穩(wěn)健地反映基因間的功能關(guān)聯(lián)
TOM = TOMsimilarity(adjacency)
# 計算 TOM 不相似度（用于聚類）
dissTOM = 1 - TOM
### (3) 基于 TOM 不相似度進行層次聚類
# 使用平均鏈接法（average linkage）對基因進行聚類
geneTree = hclust(as.dist(dissTOM), method = "average")
# 繪制初步的基因聚類樹（屏幕顯示，不保存）
plot(geneTree, xlab = "", sub = "", main = "Gene clustering on TOM-based dissimilarity",
labels = FALSE, hang = 0.04)
### (4) 使用動態(tài)樹剪枝法識別初始模塊
# cutreeDynamic 根據(jù)聚類樹的高度和內(nèi)部結(jié)構(gòu)自動切分出模塊
minModuleSize = 30 # 模塊最小基因數(shù)（與一步法的 minModuleSize 意義相同）
dynamicMods = cutreeDynamic(dendro = geneTree, # 聚類樹
distM = dissTOM, # 不相似度矩陣（輔助剪切）
deepSplit = 2, # 分裂靈敏度：0~4，值越大模塊數(shù)越多
pamRespectsDendro = FALSE, # 是否使用PAM（圍繞中心點劃分）來尊重樹結(jié)構(gòu)
minClusterSize = minModuleSize)
table(dynamicMods) # 查看每個初始模塊的基因數(shù)（數(shù)字編號）
# 將數(shù)字標(biāo)簽轉(zhuǎn)換為顏色標(biāo)簽（可視化用）
dynamicColors = labels2colors(dynamicMods)
table(dynamicColors) # 統(tǒng)計各顏色模塊基因數(shù)
# 繪制帶模塊顏色條的基因聚類樹（初始切分結(jié)果）
plotDendroAndColors(geneTree, dynamicColors, "Dynamic Tree Cut",
dendroLabels = FALSE, hang = 0.03,
addGuide = TRUE, guideHang = 0.05,
main = "Gene dendrogram and module colors")
### (5) 合并表達譜高度相似的模塊
# 計算每個模塊的“特征基因”（eigengene，即為模塊內(nèi)基因表達的第一主成分）
MEList = moduleEigengenes(datExpr, colors = dynamicColors)
MEs = MEList$eigengenes
# 計算模塊特征基因之間的不相似度（1 - Pearson相關(guān)系數(shù)）
MEDiss = 1 - cor(MEs)
# 對模塊特征基因進行層次聚類
METree = hclust(as.dist(MEDiss), method = "average")
# 設(shè)置合并閾值：若兩個模塊特征基因的相關(guān)系數(shù)達到 0.5（即不相似度 ≤ 0.5）則合并
# 此處意為：高度相似的模塊（相關(guān)系數(shù) > 0.5）將合并
MEDissThres = 0.5
# 繪制模塊特征基因聚類樹，并在合并高度處畫紅線
plot(METree, main = "Clustering of module eigengenes",
xlab = "", sub = "")
abline(h = MEDissThres, col = "red")
# 自動合并相似模塊
merge = mergeCloseModules(datExpr, dynamicColors, cutHeight = MEDissThres, verbose = 3)
# 獲取合并后的模塊顏色標(biāo)簽
mergedColors = merge$colors
# 獲取合并后的模塊特征基因
mergedMEs = merge$newMEs
# 統(tǒng)計合并后的模塊基因數(shù)
table(mergedColors)
### (6) 保存最終的可視化圖形：合并前后的模塊對比
pdf(file = "step3_stepbystep_DynamicTreeCut_genes-modules.pdf", width = 16, height = 12)
plotDendroAndColors(geneTree,
cbind(dynamicColors, mergedColors), # 左側(cè)為合并前，右側(cè)為合并后
c("Dynamic Tree Cut", "Merged dynamic"),
dendroLabels = FALSE, hang = 0.03,
addGuide = TRUE, guideHang = 0.05)
dev.off()
### 保存分布法的最終結(jié)果，供后續(xù)分析使用
# 將合并后的顏色命名為 moduleColors（與一步法輸出變量名一致）
moduleColors = mergedColors
# 生成數(shù)字標(biāo)簽（0 對應(yīng)灰色，1~N 對應(yīng)標(biāo)準顏色）
colorOrder = c("grey", standardColors(50))
moduleLabels = match(moduleColors, colorOrder) - 1
MEs = mergedMEs
save(MEs, moduleLabels, moduleColors, geneTree,
file = "step3_stepByStep_genes_modules.Rdata")
}
④ 關(guān)聯(lián)基因模塊與表型
先繪制三類圖形：模塊與表型相關(guān)性的熱圖、模塊與表型相關(guān)性的箱線圖，以及基因與模塊的相關(guān)性 vs 基因與表型的相關(guān)性的散點圖。
綜合上述結(jié)果，判斷哪些模塊與表型之間存在高度相關(guān)性，從而將基因模塊與表型聯(lián)系起來。
（以散點圖為例，當(dāng)需要查看某個離散性狀表型時，這里演示選擇的是“primed”這一表型。）
####################### 4.關(guān)聯(lián)基因模塊與表型 #####################################
# 清空環(huán)境，重新加載之前步驟生成的關(guān)鍵數(shù)據(jù)
rm(list = ls())
load(file = "step1_input.Rdata") # 包含 datExpr, datTraits, nGenes, nSamples
load(file = "step2_power_value.Rdata") # 包含 sft, power（本次未直接使用）
load(file = "step3_genes_modules.Rdata") # 包含 net, moduleColors
## 模塊與表型的相關(guān)性熱圖
if(T){
# 將分組列轉(zhuǎn)換為因子（便于后續(xù)構(gòu)建設(shè)計矩陣）
datTraits$group <- as.factor(datTraits$group)
# 創(chuàng)建設(shè)計矩陣（dummy variables），每個組別為一列，值為0/1表示樣本是否屬于該組
design <- model.matrix(~0 + datTraits$group)
colnames(design) <- levels(datTraits$group) # 列名為各個分組名稱
# 計算模塊特征基因（eigengene）：每個模塊內(nèi)基因表達值的第一主成分
MES0 <- moduleEigengenes(datExpr, moduleColors)$eigengenes
# 對特征基因進行重排序，使相似模塊的特征基因在熱圖中相鄰
MEs <- orderMEs(MES0)
# 計算模塊特征基因與表型（設(shè)計矩陣）的 Pearson 相關(guān)系數(shù)
# use = "p" 表示使用成對完全觀測（pairwise complete observations）
moduleTraitCor <- cor(MEs, design, use = "p")
# 計算相關(guān)系數(shù)的顯著性 p 值（基于 Student t 分布）
moduleTraitPvalue <- corPvalueStudent(moduleTraitCor, nSamples)
# 創(chuàng)建熱圖中的文本標(biāo)簽：顯示“相關(guān)系數(shù)\n(p值)”，保留有效數(shù)字
textMatrix <- paste0(signif(moduleTraitCor, 2), "\n(",
signif(moduleTraitPvalue, 1), ")")
dim(textMatrix) <- dim(moduleTraitCor)
# 輸出熱圖 PDF，寬度根據(jù)表型數(shù)量、高度根據(jù)模塊數(shù)量動態(tài)調(diào)整
pdf("step4_Module-trait-relationship_heatmap.pdf",
width = 2 * length(colnames(design)),
height = 0.6 * length(names(MEs)))
par(mar = c(5, 9, 3, 3)) # 設(shè)置圖形邊距（下、左、上、右）
labeledHeatmap(Matrix = moduleTraitCor,
xLabels = colnames(design), # 表型標(biāo)簽（分組名）
yLabels = names(MEs), # 模塊特征基因名
ySymbols = names(MEs), # 同上，用于圖例
colorLabels = FALSE, # 不使用顏色標(biāo)簽
colors = blueWhiteRed(50), # 顏色方案：藍-白-紅
textMatrix = textMatrix, # 顯示數(shù)值和 p 值
setStdMargins = FALSE, # 不使用默認邊距
cex.text = 0.5, # 熱圖內(nèi)文字大小
zlim = c(-1, 1), # 顏色范圍（相關(guān)系數(shù)范圍）
main = "Module-trait relationships")
dev.off()
# 保存設(shè)計矩陣，供后續(xù)分析使用
save(design, file = "step4_design.Rdata")
}
### 模塊與表型的相關(guān)性 boxplot 圖
if(T){
# 將模塊特征基因（MEs）與表型信息（datTraits）按行名（樣本名）合并
mes_group <- merge(MEs, datTraits, by = "row.names")
# 加載繪圖所需的包
library(gplots)
library(ggpubr)
library(grid)
library(gridExtra)
# 定義一個批量繪制箱線圖的函數(shù)
draw_ggboxplot <- function(data, Module = "Module", group = "group"){
ggboxplot(data, x = group, y = Module,
ylab = paste0(Module),
xlab = group,
fill = group,
palette = "jco", # 期刊配色（類似 jco 雜志風(fēng)格）
# add = "jitter", # 可添加抖動點（默認注釋掉）
legend = "") +
stat_compare_means() # 自動添加組間比較的 p 值（方法默認 t 檢驗或 Wilcoxon）
}
# 對所有模塊批量繪制箱線圖
colorNames <- names(MEs) # 模塊特征基因名稱（如 MEblue, MEbrown...）
pdf("step4_Module-trait-relationship_boxplot.pdf",
width = 7.5, height = 1.6 * ncol(MEs))
p <- lapply(colorNames, function(x) {
draw_ggboxplot(mes_group, Module = x, group = "group")
})
do.call(grid.arrange, c(p, ncol = 2)) # 將多個圖排列為每行 2 個
dev.off()
}
### 基因與模塊、表型的相關(guān)性散點圖
# 目的：檢查在某一特定性狀（如"primed"組）中，基因的模塊隸屬度（MM）與基因顯著性（GS）之間是否存在正相關(guān)。
# 提示：如果一個模塊內(nèi)的基因普遍與目標(biāo)性狀顯著相關(guān)，則該模塊很可能是與該性狀相關(guān)的關(guān)鍵模塊。
# 選擇需要分析的離散表型（需根據(jù)實際分組名稱修改）
levels(datTraits$group) # 查看所有分組名稱，選擇其中一個
choose_group <- "primed" # 示例：選擇 "primed" 作為目標(biāo)分組（此處需根據(jù)真實數(shù)據(jù)修改）
if(T){
# 提取模塊名稱（去掉特征基因名中的 "ME" 前綴）
modNames <- substring(names(MEs), 3) # 例如 "MEblue" -> "blue"
### 計算模塊與基因的相關(guān)性矩陣 —— 模塊成員關(guān)系（Module Membership, MM）
# MM 值為每個基因的表達譜與對應(yīng)模塊特征基因的 Pearson 相關(guān)系數(shù)
geneModuleMembership <- as.data.frame(cor(datExpr, MEs, use = "p"))
# 計算 MM 的 p 值
MMPvalue <- as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples))
# 修改列名，使其包含模塊顏色（例如 MMblue）
names(geneModuleMembership) <- paste0("MM", modNames)
names(MMPvalue) <- paste0("p.MM", modNames)
### 計算性狀與基因的相關(guān)性矩陣 —— 基因顯著性（Gene Significance, GS）
# GS 為每個基因的表達值與目標(biāo)表型（0/1 啞變量）的 Pearson 相關(guān)系數(shù)
## 對于非連續(xù)型性狀，應(yīng)將其轉(zhuǎn)換為 0-1 矩陣，這里直接使用先前創(chuàng)建設(shè)計矩陣中的某一列（代表是否屬于 choose_group）
trait <- as.data.frame(design[, choose_group])
geneTraitSignificance <- as.data.frame(cor(datExpr, trait, use = "p"))
GSPvalue <- as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples))
names(geneTraitSignificance) <- "GS"
names(GSPvalue) <- "GS"
### 可視化：對每個模塊，繪制其內(nèi)部基因的 |MM| 與 |GS| 的散點圖
# 可以選擇特定模塊，也可以遍歷全部模塊
# selectModule <- c("blue","green","purple","grey") # 手動指定關(guān)注的模塊
selectModule <- modNames # 預(yù)設(shè)為分析全部模塊
pdf("step4_gene-Module-trait-significance.pdf", width = 7, height = 1.5 * ncol(MEs))
par(mfrow = c(ceiling(length(selectModule) / 2), 2)) # 多圖排布：每行2個
for(module in selectModule){
column <- match(module, selectModule) # 獲取該模塊在 selectModule 中的索引位置
print(module) # 輸出當(dāng)前模塊名，便于控制臺查看進度
moduleGenes <- (moduleColors == module) # 邏輯向量：屬于該模塊的基因
# 繪制 |MM| 與 |GS| 的散點圖，并添加擬合線（由 verboseScatterplot 自動添加）
verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
abs(geneTraitSignificance[moduleGenes, 1]),
xlab = paste("Module Membership in", module, "module"),
ylab = "Gene significance for trait",
main = paste("Module membership vs. gene significance\n"),
cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2,
col = module) # 點的顏色與模塊顏色一致
}
dev.off()
}
⑤ WGCNA可視化：TOMplot 、Eigengene-adjacency-heatmap
在 WGCNA 中，標(biāo)配的熱圖 TOMplot（或稱網(wǎng)絡(luò)熱圖）展示的是全基因范圍內(nèi)的拓撲重疊矩陣（TOM），顏色越深代表基因間的共表達相關(guān)性越強。
另一張熱圖——特征基因鄰接熱圖（Eigengene-adjacency heatmap）則用于呈現(xiàn)各基因模塊之間的相關(guān)性，同時還可以把表型數(shù)據(jù)添加進去，從而直觀地看出哪些模塊與特定表型相關(guān)性較高。
######################### 5. WGCNA可視化：TOMplot Eigengene-adjacency-heatmap ##################################
# 清空環(huán)境并加載之前分析生成的關(guān)鍵數(shù)據(jù)
rm(list = ls())
load(file = 'step1_input.Rdata') # 包含 datExpr, datTraits, nGenes, nSamples
load(file = "step2_power_value.Rdata") # 包含 sft, power（軟閾值冪次）
load(file = "step3_genes_modules.Rdata") # 包含 net, moduleColors
load(file = "step4_design.Rdata") # 包含 design（分組設(shè)計矩陣）
# ############## 第一部分：基于TOM矩陣的全網(wǎng)絡(luò)熱圖 ##############
if(T){
# 從表達數(shù)據(jù)直接計算拓撲重疊矩陣（TOM）
# TOMsimilarityFromExpr 內(nèi)部會自動執(zhí)行：相關(guān)性計算 → 軟閾值冪次轉(zhuǎn)換 → TOM計算
# 參數(shù) power 即為之前確定的軟閾值
TOM = TOMsimilarityFromExpr(datExpr, power = power)
# 計算 TOM 不相似度（用于距離度量）
dissTOM = 1 - TOM
## 繪制所有基因的 TOM 熱圖（通常用于展示網(wǎng)絡(luò)整體結(jié)構(gòu)，但基因較多時文件很大）
if(T){
# 從 net 對象中提取基因聚類樹（一步法的結(jié)果）
geneTree = net$dendrograms[[1]]
# 對 TOM 矩陣進行冪次轉(zhuǎn)換（此處為7次方），旨在增強對比度，使高連接強度的區(qū)域更明顯
# 冪次越大，弱連接趨近于0，強連接保留較高值，熱圖顏色分界更清晰
plotTOM = dissTOM^7
diag(plotTOM) = NA # 將對角線（基因自身連接）設(shè)為缺失值，避免干擾視覺
# 輸出為 PNG 格式（適合屏幕查看，但可能不如 PDF 清晰）
png("step5_TOMplot_Network-heatmap.png", width = 800, height = 600)
# TOMplot 將 TOM 熱圖、基因聚類樹和模塊顏色條整合在一張圖中
TOMplot(plotTOM, geneTree, moduleColors,
col = gplots::colorpanel(250, 'red', "orange", 'lemonchiffon'),
main = "Network heatmap plot")
dev.off()
}
### 為了節(jié)省繪圖時間和文件大小，可以隨機挑選部分基因進行測試（當(dāng)前未執(zhí)行）
if(F){
nSelect = 0.1 * nGenes # 選擇10%的基因
set.seed(123) # 設(shè)置隨機種子，保證結(jié)果可重復(fù)
select = sample(nGenes, size = nSelect)
selectTOM = dissTOM[select, select] # 子矩陣
selectTree = hclust(as.dist(selectTOM), method = "average")
selectColors = moduleColors[select]
plotDiss = selectTOM^7
diag(plotDiss) = NA
pdf("step5_select_TOMplot_Network-heatmap.pdf", width = 8, height = 6)
TOMplot(plotDiss, selectTree, selectColors,
col = gplots::colorpanel(250, 'red', "orange", 'lemonchiffon'),
main = "Network heatmap plot of selected gene")
dev.off()
}
}
# ############## 第二部分：模塊特征基因（eigengene）相關(guān)性網(wǎng)絡(luò) ##############
# 這部分展示模塊之間的關(guān)系，以及模塊與表型之間的關(guān)系
if(T){
# 計算每個模塊的特征基因（eigengene，即模塊內(nèi)表達譜的第一主成分）
MEs = moduleEigengenes(datExpr, moduleColors)$eigengenes
# 對特征基因進行排序，使相似的模塊在后續(xù)熱圖中相鄰
MET = orderMEs(MEs)
# 可選擇將表型數(shù)據(jù)（連續(xù)或分類）添加到模塊特征基因集中，一起展示相關(guān)性
if(T){
## 對于連續(xù)型性狀，可直接添加數(shù)值列，例如：MET = orderMEs(cbind(MEs, datTraits$groupNo))
## 對于離散性狀，需要使用啞變量（已經(jīng)保存在 design 中）
# 此處示例將 design 中的第三列（假設(shè)代表 "primed" 組）合并進去
# 注意 design 的列順序取決于分組水平的順序，需根據(jù)實際情況調(diào)整索引
design
primed = as.data.frame(design[, 3]) # 提取第三列（例如 "primed" 組）
names(primed) = "primed"
# 將表型列與模塊特征基因合并，并重新排序
MET = orderMEs(cbind(MEs, primed))
}
# 繪制模塊特征基因的聚類樹和熱圖（展示模塊間相關(guān)性及模塊與表型的相關(guān)性）
pdf("step5_module_cor_Eigengene-dendrogram.pdf", width = 8, height = 10)
plotEigengeneNetworks(MET,
setLabels = "", # 圖中標(biāo)題留空
marDendro = c(0, 4, 1, 4), # 聚類樹部分的邊距（下、右、上、左，單位：行文本高度）
marHeatmap = c(5, 5, 1, 2), # 熱圖部分的邊距
cex.lab = 0.8, # 坐標(biāo)軸標(biāo)簽字號
xLabelsAngle = 90) # X軸標(biāo)簽旋轉(zhuǎn)角度（90度垂直顯示）
dev.off()
}
從上圖可以發(fā)現(xiàn)，primed表型與"brown"模塊聚類是最接近的
⑥ 對感興趣模塊的基因進行GO分析
這里選取了"black","pink","brown"模塊進行GO分析
#################### 6. 選擇感興趣基因模塊進行GO分析 ####################
# 清空環(huán)境，加載前面步驟生成的必需數(shù)據(jù)
rm(list = ls())
load(file = 'step1_input.Rdata') # 包含 datExpr（表達矩陣）等
load(file = "step2_power_value.Rdata") # 包含 power 值
load(file = "step3_genes_modules.Rdata") # 包含 net, moduleColors
load(file = "step4_design.Rdata") # 包含 design（分組設(shè)計矩陣）
### 條件設(shè)置（需根據(jù)實際物種修改）
OrgDb = "org.Mm.eg.db" # 物種注釋包：小鼠為 org.Mm.eg.db，人類為 org.Hs.eg.db
genetype = "SYMBOL" # 基因名類型："SYMBOL" 或 "ENSEMBL"（需與表達矩陣行名一致）
table(moduleColors) # 查看各個模塊包含的基因數(shù)量，幫助選擇感興趣的模塊
choose_module <- c("black", "pink", "brown") # 手動選擇要分析的模塊（按顏色名）
if(T){
# 加載 GO 分析相關(guān)的 R 包
library(clusterProfiler) # 富集分析主包
library(org.Mm.eg.db) # 小鼠基因注釋（若改為人類則需加載 org.Hs.eg.db）
library(org.Hs.eg.db) # 人類基因注釋（此處雙加載，實際按需選擇一個或注釋掉另一個）
### 1. 建立基因與模塊的對應(yīng)關(guān)系表
gene_module <- data.frame(gene = colnames(datExpr), # 基因名（來自表達矩陣列名）
module = moduleColors) # 對應(yīng)模塊顏色
# 保存為 CSV 文件，便于后續(xù)查閱
write.csv(gene_module, file = "step6_gene_moduleColors.csv", row.names = FALSE, quote = FALSE)
### 2. 將基因名轉(zhuǎn)換為 ENTREZ ID（富集分析通常需要 ENTREZ ID）
# bitr 函數(shù)用于 ID 轉(zhuǎn)換：從 genetype（如 SYMBOL）轉(zhuǎn)換為 "ENTREZID"
tmp <- bitr(gene_module$gene,
fromType = genetype, # 輸入 ID 類型
toType = "ENTREZID", # 輸出 ID 類型
OrgDb = OrgDb) # 對應(yīng)的物種注釋包
# 將轉(zhuǎn)換后的 ENTREZ ID 與原模塊信息合并
gene_module_entrz <- merge(tmp, gene_module, by.x = genetype, by.y = "gene")
### 3. 篩選出指定模塊的基因（用于后續(xù)富集分析）
choose_gene_module_entrz <- gene_module_entrz[gene_module_entrz$module %in% choose_module, ]
### 4. 運行 GO 富集分析（按模塊分組進行比較）
# compareCluster 可以對多個基因列表（按模塊分組）同時進行 GO 富集分析
formula_res <- compareCluster(
ENTREZID ~ module, # 公式：基因列 ~ 分組列（模塊名）
data = choose_gene_module_entrz, # 包含 ENTREZID 和 module 列的數(shù)據(jù)框
fun = "enrichGO", # 使用的富集函數(shù)（GO 分析）
OrgDb = OrgDb, # 注釋數(shù)據(jù)庫
ont = "BP", # GO 子本體：BP（生物學(xué)過程），也可選 MF、CC 或 "ALL"
pAdjustMethod = "BH", # 多重檢驗校正方法（Benjamini-Hochberg）
pvalueCutoff = 0.25, # p 值閾值（較寬松，可根據(jù)需要調(diào)整）
qvalueCutoff = 0.25 # q 值閾值（FDR 校正后）
)
### 5. 精簡富集結(jié)果，去除冗余的 GO 條目
# simplify 函數(shù)可根據(jù)語義相似度合并相似的 GO 術(shù)語
lineage1_ego <- simplify(
formula_res,
cutoff = 0.5, # 相似度閾值（0.5 表示中等嚴格）
by = "p.adjust", # 選擇保留 p.adjust 最小的條目
select_fun = min # 當(dāng)合并時，選擇最小 p.adjust 的條目為代表
)
# 保存富集分析的中間結(jié)果和最終結(jié)果
save(gene_module, formula_res, lineage1_ego, file = "step6_module_GO_term.Rdata")
# 將精簡后的結(jié)果寫入 CSV 文件（便于在 Excel 中查看）
write.csv(lineage1_ego@compareClusterResult,
file = "step6_module_GO_term.csv")
### 6. 繪制點圖（dotplot）展示各模塊顯著富集的 GO 條目
dotp <- dotplot(lineage1_ego,
showCategory = 10, # 每個模塊展示最多 10 個 GO 條目
includeAll = TRUE, # 包含有重疊（但未顯著）的結(jié)果（通常保持默認）
label_format = 90) # 標(biāo)簽過長時自動換行的字符寬度
# 保存為 PDF 文件
ggsave(dotp, filename = "step6_module_GO_term.pdf",
width = 12, height = 15)
}
⑦ 感興趣基因模塊繪制熱圖
這里選取了與primed表型相關(guān)度很高的brown模塊繪制熱圖
############################### 7.感興趣基因模塊繪制熱圖 ######################################
# 清空環(huán)境并加載所需數(shù)據(jù)
rm(list = ls())
load(file = 'step1_input.Rdata') # 包含 datExpr（表達矩陣）、datTraits（樣本表型）
load(file = "step3_genes_modules.Rdata") # 包含 net, moduleColors（基因所屬模塊顏色）
# 查看各模塊基因數(shù)量，確認要繪制的模塊名稱
table(moduleColors)
# 指定感興趣的模塊（需與 table 輸出中的顏色名一致）
module = "brown"
### 對感興趣模塊內(nèi)的基因繪制表達熱圖
if(T){
# 提取屬于該模塊的所有基因的表達數(shù)據(jù)（datExpr 行為樣本，列為基因）
dat = datExpr[, moduleColors == module]
# 加載 pheatmap 包（如果沒有安裝，需要先 install.packages("pheatmap")）
library(pheatmap)
# 對每個基因（列）進行 Z-score 標(biāo)準化（scale），然后轉(zhuǎn)置為行為基因、列為樣本的矩陣
# scale 默認對列進行操作，dat 的列是基因，所以標(biāo)準化后每列均值為0，標(biāo)準差為1
# t() 轉(zhuǎn)置后：行 = 基因，列 = 樣本
n = t(scale(dat)) # 標(biāo)準化后矩陣
# 可選操作：將表達值截斷到 [-2, 2] 范圍，避免個別極值影響顏色對比（當(dāng)前被注釋）
# n[n > 2] = 2
# n[n < -2] = -2
# 獲取樣本分組信息（用于熱圖上方的樣本注釋條）
group_list = datTraits$group
# 構(gòu)建注釋數(shù)據(jù)框：行名與轉(zhuǎn)置后矩陣的列名（樣本名）對應(yīng)
ac = data.frame(g = group_list)
rownames(ac) = colnames(n)
# 使用 pheatmap 繪制熱圖并保存為 PDF
pheatmap::pheatmap(n,
fontsize = 8, # 全局字體大小
show_colnames = TRUE, # 顯示列名（樣本名）
show_rownames = FALSE, # 不顯示行名（基因名太多，避免擁擠）
cluster_cols = TRUE, # 對樣本進行列聚類
annotation_col = ac, # 添加樣本注釋條（分組信息）
width = 8, # 輸出圖片寬度（英寸）
height = 6, # 輸出圖片高度（英寸）
angle_col = 45, # 列名（樣本名）傾斜角度
main = paste0("module_", module, "-gene heatmap"), # 主標(biāo)題
filename = paste0("step7_module_", module, "_Gene-heatmap.pdf"))
}
⑧ 導(dǎo)出基因至 VisANT 或 cytoscape
將感興趣模塊 "brown"的基因?qū)С?VisANT或cytoscape
################### 8.感興趣模塊基因?qū)С?VisANT or cytoscape ######################
# 清空環(huán)境并加載分析所需的數(shù)據(jù)
rm(list = ls())
load(file = 'step1_input.Rdata') # 包含 datExpr（表達矩陣）、datTraits 等
load(file = "step2_power_value.Rdata") # 包含 power（軟閾值冪次）
load(file = "step3_genes_modules.Rdata") # 包含 net, moduleColors
# 指定要導(dǎo)出的模塊顏色名稱（需與 moduleColors 中的顏色對應(yīng)）
module = "brown"
if(T){
### 1. 提取感興趣模塊的基因名
gene <- colnames(datExpr) # 所有基因名
inModule <- moduleColors == module # 邏輯向量：TRUE 表示屬于該模塊
modgene <- gene[inModule] # 該模塊內(nèi)的基因名列表
### 2. 計算該模塊對應(yīng)的基因關(guān)系矩陣（拓撲重疊矩陣 TOM）
# TOMsimilarityFromExpr 基于表達數(shù)據(jù)和軟閾值重新計算整個網(wǎng)絡(luò)的 TOM
# （若數(shù)據(jù)量很大，可考慮提前保存整個 TOM，避免重復(fù)計算）
TOM <- TOMsimilarityFromExpr(datExpr, power = power)
# 提取子矩陣：僅包含該模塊基因之間的 TOM 值
modTOM <- TOM[inModule, inModule]
# 設(shè)置行列名稱（為基因名，便于導(dǎo)出后識別）
dimnames(modTOM) <- list(modgene, modgene)
### 3. 篩選連接度最大的前 N 個基因（例如 Top 100），以減少網(wǎng)絡(luò)復(fù)雜度
nTop = 100
# softConnectivity 計算每個基因在給定子集中的連接度（即與其他基因的 TOM 值之和）
IMConn = softConnectivity(datExpr[, modgene])
# 選取連接度排名前 nTop 的基因（rank(-IMConn) 表示從大到小排序，取前 nTop 個）
top = (rank(-IMConn) <= nTop)
# 進一步篩選 TOM 子矩陣，只保留這些高連接度的基因
filter_modTOM <- modTOM[top, top]
# 4. 導(dǎo)出為 VisANT 輸入格式
# VisANT 是一個網(wǎng)絡(luò)可視化軟件，exportNetworkToVisANT 可生成包含邊信息的文本文件
vis <- exportNetworkToVisANT(filter_modTOM,
file = paste("step8_visANTinput-", module, ".txt", sep = ""),
weighted = TRUE, # 輸出權(quán)重信息（即 TOM 值）
threshold = 0) # 閾值設(shè)為 0 表示保留所有邊
# 注意：VisANT 格式通常只需要邊列表，節(jié)點信息可自動從邊列表提取
# 5. 導(dǎo)出為 Cytoscape 輸入格式（更主流）
# Cytoscape 需要的文件：邊文件（edgeFile）和節(jié)點文件（nodeFile）
cyt <- exportNetworkToCytoscape(filter_modTOM,
edgeFile = paste("step8_CytoscapeInput-edges-", paste(module, collapse="-"), ".txt", sep=""),
nodeFile = paste("step8_CytoscapeInput-nodes-", paste(module, collapse="-"), ".txt", sep=""),
weighted = TRUE, # 輸出權(quán)重（后可用作邊粗細）
threshold = 0.15, # 權(quán)重閾值：只保留 TOM > 0.15 的邊（可調(diào)整以簡化網(wǎng)絡(luò)）
nodeNames = modgene[top], # 節(jié)點名稱（高連接度的基因）
nodeAttr = moduleColors[inModule][top]) # 節(jié)點屬性（此處為模塊顏色，實際所有節(jié)點屬于同一模塊，顏色相同）
}
