加權(quán)基因共表達網(wǎng)絡(luò)分析（WGCNA, Weighted gene co-expression network analysis）是一種用來描述不同樣品之間基因關(guān)聯(lián)模式的系統(tǒng)生物學方法争占，可以用來鑒定高度協(xié)同變化的基因集, 并根據(jù)基因集的內(nèi)連性和基因集與表型之間的關(guān)聯(lián)鑒定候補生物標記基因或治療靶點输硝。相比于只關(guān)注差異表達的基因痒钝，WGCNA利用數(shù)千或近萬個變化最大的基因或全部基因的信息識別感興趣的基因集，并與表型進行顯著性關(guān)聯(lián)分析锤灿。一是充分利用了信息钝鸽，二是把數(shù)千個基因與表型的關(guān)聯(lián)轉(zhuǎn)換為數(shù)個基因集與表型的關(guān)聯(lián)嘉抓，免去了多重假設(shè)檢驗校正的問題毡泻。
專有名詞：
? 共表達網(wǎng)絡(luò)：定義為加權(quán)基因網(wǎng)絡(luò)。點代表基因媳瞪，邊代表基因表達相關(guān)性
? Module(模塊)：高度內(nèi)連的基因集骗炉。在無向網(wǎng)絡(luò)中，模塊內(nèi)是高度相關(guān)的基因蛇受；在有向網(wǎng)絡(luò)中句葵，模塊內(nèi)是高度正相關(guān)的基因。把基因聚類成模塊后龙巨，可以對每個模塊進行三個層次的分析：1. 功能富集分析查看其功能特征是否與研究目的相符笼呆；2. 模塊與性狀進行關(guān)聯(lián)分析，找出與關(guān)注性狀相關(guān)度最高的模塊旨别；3. 模塊與樣本進行關(guān)聯(lián)分析诗赌，找到樣品特異高表達的模塊
? 連接度(Connectivity)：類似于網(wǎng)絡(luò)中 "度"(degree)的概念。每個基因的連接度是與其相連的基因的邊屬性之和
? 鄰接矩陣(Adjacency matrix)：基因和基因之間的加權(quán)相關(guān)性值構(gòu)成的矩陣
? Intramodular connectivity：模內(nèi)連通性測量秸弛，給定基因相對于特定模塊的基因如何連接或共表達铭若，判斷基因所屬關(guān)系
? Module eigengene E：給定模型的第一主成分洪碳，代表整個模型的基因表達譜，可很好的用一個向量代替了一個矩陣叼屠，方便后期計算
? Eigengene signi?cance：當獲得微陣列樣品特征y（例如病例對照狀態(tài)或體重）時瞳腌，可以將模塊特征基因與此結(jié)果相關(guān)聯(lián)， 相關(guān)系數(shù)稱為特征基因顯著性
? Module Membership KME ：對于每個基因镜雨，通過將其基因表達譜與給定模塊的模塊本征基因相關(guān)聯(lián)來定義模塊成員的“模糊”度量嫂侍。如，用來測量基因
與藍色模塊特征基因的相關(guān)程度荚坞。若
接近0挑宠，則第
? Gene signi?cance GS：的絕對值越高诡挂，第
? Module signi?cance：給定模塊中所有基因的平均絕對基因顯著性的度量。 當將基因顯著性定義為基因表達與外部性狀y的相關(guān)性時璃俗，此度量往往與模塊特征基因與y的相關(guān)性高度相關(guān)
? TOM (Topological overlap matrix)：把鄰接矩陣轉(zhuǎn)換為拓撲重疊矩陣奴璃，以降低噪音和假相關(guān)，獲得的新距離矩陣城豁，這個信息可拿來構(gòu)建網(wǎng)絡(luò)或繪制TOM圖
? Hub gene：關(guān)鍵基因 (連接度最多或連接多個模塊的基因)
? 軟閾值：相關(guān)性值進行冪次運算冪次的值也就是軟閾值;
WGCNA-1
1. 數(shù)據(jù)輸入溺健、清洗和預(yù)處理
1.1 加載基因表達數(shù)據(jù)
library(WGCNA)
dim(FPKM) # 仍選取TCGA- LUAD的FPKM數(shù)據(jù)
## [1] 59427 600
## 篩選中位絕對偏差前75%的基因，至少MAD大于0.01
## 篩選后會降低運算量钮蛛，也會失去部分信息
## 也可不做篩選，使MAD大于0即可
m.mad <- apply(FPKM, 1, mad)
dataExprVar <- FPKM[which(m.mad > max(quantile(m.mad, probs = seq(0, 1, 0.25))[2], 0.01)),]
# 過濾基因剖膳，選取絕對中位差Top5000的基因
# dataExprVar <- t(exp[order(apply(FPKM, 1, mad), decreasing = T)[1:5000],])
dim(dataExprVar)
## [1] 30976 600
## 轉(zhuǎn)換為樣品在行魏颓，基因在列的矩陣
LUAD_Expr0 <- as.data.frame(t(dataExprVar))
★ 注意：并不推薦使用差異基因作為輸入矩陣，通過差異表達基因過濾將會導致一個（或幾個高度相關(guān)的）基因聚成一個模塊吱晒，同時甸饱，也破壞了無標度拓撲的假設(shè)，所以通過無標度拓撲擬合來選擇軟閾值的將會失敗
1.2 數(shù)據(jù)清洗
## 檢測缺失值
gsg = goodSamplesGenes(LUAD_Expr0, verbose = 3)
gsg$allOK;
## 1] TRUE
# 若樣本檢查語句返回TRUE仑濒，那么沒有缺失值叹话。否則，可通過以下代碼來移除缺失值
if (!gsg$allOK){
# Optionally, print the gene and sample names that were removed:
if (sum(!gsg$goodGenes)>0)
printFlush(paste("Removing genes:",
paste(names(LUAD_Expr0)[!gsg$goodGenes], collapse = ",")));
if (sum(!gsg$goodSamples)>0)
printFlush(paste("Removing samples:",
paste(rownames(LUAD_Expr0)[!gsg$goodSamples], collapse = ",")));
# Remove the offending genes and samples from the data:
LUAD_Expr0 = LUAD_Expr0[gsg$goodSamples, gsg$goodGenes]
}
## 對樣本進行聚類（與隨后的基因聚類相比）墩瞳，看看是否有明顯的異常值
# hclusts聚類算法, dist計算基因之間的距離
sampleTree = hclust(dist(LUAD_Expr0), method = "average");
# 設(shè)置文字大小
par(cex = 0.2);
# 設(shè)置圖像邊距c(bottom, left, top, right)
par(mar = c(0,4,2,0))
# 畫圖 main標題驼壶，sub子標題，xlab x軸標題喉酌，cex.lab標題字體大小热凹，cex.axis坐標軸刻度大小泵喘，cex.main主標題字體
plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5, cex.axis = 1.5, cex.main = 5)
# 在上圖上畫紅線
abline(h = 1e+05, col = "red");
# Determine cluster under the line
# 剪枝算法，cutHeight 修剪樹枝的高度 minSize集群最小數(shù)
clust = cutreeStatic(sampleTree, cutHeight = 1e+05, minSize = 5)
# 剪枝結(jié)果
table(clust)
## clust
## 0 1 2
## 2 591 7
# clust 1 contains the samples we want to keep
keepSamples = (clust==1)
# 符合要求的數(shù)據(jù)
LUAD_Expr = LUAD_Expr0[keepSamples, ]
# 提取行和列
nSamples = nrow(LUAD_Expr)
nGenes = ncol(LUAD_Expr)
WGCNA-2
★ 此處去除左側(cè)9個樣本
1.3 加載臨床特征數(shù)據(jù)
# 讀取臨床數(shù)據(jù)
LUAD_clin0 <- jsonlite::fromJSON("/data/shumin/Cibersort/LUAD/clinical.cart.2023-06-08.json")
n = length(LUAD_clin0$diagnoses)
# 提取數(shù)據(jù)
case_id = classfication_of_tumor = c(rep(0, n))
tumor_stage = gender = c(rep(0, n))
year_to_birth = year_to_death = c(rep(0, n))
year_to_diagnosis = days_to_death = c(rep(0, n))
age = deadORlive = race = alcohol = smoked = c(rep(0, n))
for (i in 1:n) {
case_id[i] = LUAD_clin0$case_id[[I]]
classfication_of_tumor[i] = LUAD_clin0$diagnoses[[i]]$classification_of_tumor
tumor_stage[i] = LUAD_clin0$diagnoses[[i]]$tumor_grade
gender[i] = LUAD_clin0$demographic$gender[[I]]
year_to_birth[i] = ifelse(
is.null(LUAD_clin0$demographic$year_of_birth[[i]]),
"notReport",
LUAD_clin0$demographic$year_of_birth[[I]]
)
year_to_death[i] = ifelse(
is.null(LUAD_clin0$demographic$year_of_death[[i]]),
"notReport",
LUAD_clin0$demographic$year_of_death[[I]]
)
year_to_diagnosis[i] = ifelse(
is.null(LUAD_clin0$diagnoses[[i]]$year_of_diagnosis),
"notReport",
LUAD_clin0$diagnoses[[i]]$year_of_diagnosis
)
days_to_death[i] = ifelse(
is.null(LUAD_clin0$demographic$days_to_death[[i]]),
"notReport",
LUAD_clin0$demographic$days_to_death[[I]]
)
age[i] = ifelse(
is.null(LUAD_clin0$demographic$age_at_index[[i]]),
"notReport",
LUAD_clin0$demographic$age_at_index[[I]]
)
deadORlive[i] = ifelse(
is.null(LUAD_clin0$demographic$vital_status[[i]]),
"notReport",
LUAD_clin0$demographic$vital_status[[I]]
)
race[i] = ifelse(
is.null(LUAD_clin0$demographic$race[[i]]),
"notReport",
LUAD_clin0$demographic$race[[I]]
)
alcohol[i] = ifelse(
is.null(LUAD_clin0$exposures[[i]]$alcohol_history),
"notReprot",
LUAD_clin0$exposures[[i]]$alcohol_history
)
smoked[i] = ifelse(
is.null(LUAD_clin0$exposures[[i]]$years_smoked),
"notReport",
LUAD_clin0$exposures[[i]]$years_smoked
)
}
LUAD_clin1 <- data.frame(
case_id,
classfication_of_tumor,
tumor_stage,
gender,
year_to_birth,
year_to_death,
year_to_diagnosis,
days_to_death,
age,
deadORlive,
race,
alcohol,
smoked
)
head(LUAD_clin1)
## case_id classfication_of_tumor tumor_stage gender year_to_birth year_to_death year_to_diagnosis days_to_death age deadORlive race alcohol smoked
## 1 a3fd20b2-e001-44ab-9716-754e5ae70808 not reported not reported female 1935 NA 2012 NA 77 Alive white Not Reported NA
## 2 25d4ea9e-f773-4f11-bac9-64efdad73211 not reported not reported female 1950 NA 2009 NA 59 Alive white Not Reported 45
## 3 a52e99d6-a61a-439d-b0b1-ca7a0eabcb04 not reported not reported female 1969 NA 2011 NA 42 Alive white Not Reported 2
## 4 27fceec1-3298-4cdd-a4e6-8f5cf34604f0 not reported not reported male 1957 NA 2010 NA 53 Alive black or african american Not Reported NA
## 5 2923e404-38f2-437a-b57e-23401fbe0273 not reported not reported male 1965 NA 2011 NA 46 Alive white Not Reported NA
## 6 a65700c2-e58c-4fd4-aeb1-5686b8f4d212 not reported not reported male 1935 NA 2009 NA 74 Alive white Not Reported NA
# 讀取meta數(shù)據(jù)（用于id轉(zhuǎn)換）
LUAD_meta <- jsonlite::fromJSON("/data/shumin/Cibersort/LUAD/metadata.cart.2023-06-08.json")
sample_id <- sapply(LUAD_meta$associated_entities, function(x){x[,1]})
case_id <- sapply(LUAD_meta$associated_entities, function(x){x[,3]})
anno_data <- data.frame(sample_id, case_id)
head(anno_data)
## sample_id case_id
## 1 TCGA-44-8120-01A-11R-2241-07 83e38dbd-edab-47f2-b19f-6ea38fc6bece
## 2 TCGA-99-8025-01A-11R-2241-07 84c3ba70-afa7-4b69-be69-7ec8d6022c56
## 3 TCGA-50-6594-01A-11R-1755-07 8504fd86-a70a-4cba-9ec8-25c9e60ca549
## 4 TCGA-78-7161-01A-11R-2039-07 86faf16c-56fd-4b7b-a6b2-b4c83bec93e1
## 5 TCGA-55-7903-01A-11R-2170-07 77c4dbb2-eceb-4e0d-bcde-63dc817d5f35
## 6 TCGA-38-4632-01A-01R-1755-07 875333ab-9048-462d-aaa2-693ad127e3cc
LUAD_clin2 <- merge(anno_data, LUAD_clin1, by = "case_id")
LUAD_clin3 <- LUAD_clin2[, c("sample_id","gender","year_to_birth","year_to_death","year_to_diagnosis","days_to_death", "age","deadORlive","race","smoked")] # 選取后續(xù)分析需要的列
LUAD_clin3$gender <- ifelse(LUAD_clin3$gender == 'female', 1, 2)
LUAD_clin3$deadORlive <- ifelse(LUAD_clin3$deadORlive == 'Alive', 0, 1)
LUAD_clin3$race <- ifelse(LUAD_clin3$race == 'not reported', 0,
ifelse (LUAD_clin3$race == 'american indian or alaska native', 1,
ifelse(LUAD_clin3$race == 'asian', 2,
ifelse(LUAD_clin3$race == 'black or african american', 3, 4))))
# 形成一個類似于表達數(shù)據(jù)的數(shù)據(jù)框架般妙，以保存臨床特征
# 提取行名
LUAD_Samples = rownames(LUAD_Expr)
# 數(shù)據(jù)匹配 返回匹配行
Trait_Rows = match(LUAD_Samples, LUAD_clin3$sample_id);
# 提取指定要求行
data_Traits = LUAD_clin3[Trait_Rows, -1];
# 提取行名
rownames(data_Traits) = LUAD_clin3[Trait_Rows, 1];
# 垃圾回收
collectGarbage();
# Re-cluster samples
# 畫聚類圖
LUAD_Tree = hclust(dist(LUAD_Expr), method = "average")
# Convert traits to a color representation: white means low, red means high, grey means missing entry
# 畫表型的熱圖
# 將特征轉(zhuǎn)換為顏色表示：白色表示低纪铺，紅色表示高，灰色表示缺少條目
# 如果signed為true 以綠色開頭代表最大負值碟渺，以白色開頭代表零附近的值鲜锚，然后變?yōu)榧t色代表正值
traitColors = numbers2colors(data_Traits, signed = FALSE);
# Plot the sample dendrogram and the colors underneath.
# 繪制出樹狀圖和下面的顏色
plotDendroAndColors(LUAD_Tree, traitColors, groupLabels = names(data_Traits), dendroLabels = FALSE, main = "LUAD dendrogram and trait heatmap")
WGCNA-3
2. 建設(shè)表達網(wǎng)絡(luò)與模塊檢測
此步驟是使用WGCNA方法進行所有網(wǎng)絡(luò)分析的基礎(chǔ)。WGCNA提出三種不同的方法構(gòu)建網(wǎng)絡(luò)并識別模塊：
? 使用方便的一步網(wǎng)絡(luò)結(jié)構(gòu)和模塊檢測功能苫拍，適合希望以最小努力達到結(jié)果的用戶芜繁；
? 為希望使用定制/替代方法進行實驗的用戶逐步構(gòu)建網(wǎng)絡(luò)和模塊檢測；
? 一種自動分塊網(wǎng)絡(luò)結(jié)構(gòu)和模塊檢測方法怯疤，適用于希望分析太大而無法同時分析的數(shù)據(jù)集的用戶浆洗。
2.1 自動一步構(gòu)建網(wǎng)絡(luò)與模塊檢測
2.1.1 軟閾值的選擇：網(wǎng)絡(luò)拓撲分析
構(gòu)建一個加權(quán)基因網(wǎng)絡(luò)需要選擇軟閾值冪β來計算鄰接矩陣權(quán)重參數(shù)，即將基因間的相關(guān)系數(shù)進行乘方運算來表征其相關(guān)性集峦，首先需要確定乘方的值
# Choose a set of soft-thresholding powers
# 給出候選的β值伏社，c(1:10)表示1到10；seq(from = 12, to = 20, by = 2)表示從12開始間隔兩個數(shù)到20
powers = c(c(1:10), seq(from = 12, to = 20, by = 2))
powers
# Call the network topology analysis function 調(diào)用網(wǎng)絡(luò)拓撲分析函數(shù)
# verbose表示輸出結(jié)果詳細程度
sft = pickSoftThreshold(LUAD_Expr, powerVector = powers, verbose = 0);
# sft這中保存了每個powers值計算出來的網(wǎng)絡(luò)特征,其中powerEstimate就是最佳power值塔淤，fitIndices保存了每個power對應(yīng)的網(wǎng)絡(luò)的特征摘昌。
str(sft)
## List of 2
## $ powerEstimate: num 2
## $ fitIndices :'data.frame': 15 obs. of 7 variables:
## ..$ Power : num [1:15] 1 2 3 4 5 6 7 8 9 10 ...
## ..$ SFT.R.sq : num [1:15] 0.123 0.916 0.916 0.932 0.938 ...
## ..$ slope : num [1:15] -0.295 -1 -0.996 -0.964 -0.954 ...
## ..$ truncated.R.sq: num [1:15] 0.393 0.921 0.94 0.933 0.925 ...
## ..$ mean.k. : num [1:15] 4171 1288 613 364 244 ...
## ..$ median.k. : num [1:15] 3579.7 676.1 173.7 54 19.4 ...
## ..$ max.k. : num [1:15] 9260 5325 3783 2938 2403 ...
# Plot the results 結(jié)果繪圖
# 設(shè)置圖的顯示一行兩列
par(mfrow = c(1,2));
cex1 = 0.9;
# Scale-free topology fit index as a function of the soft-thresholding power
# 生成閾值和網(wǎng)絡(luò)的特征之間的關(guān)系函數(shù)
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");
# this line corresponds to using an R^2 cut-off of h
abline(h = 0.90, col = "red")
# sft$fitIndices 保存了每個power構(gòu)建的相關(guān)性網(wǎng)絡(luò)中的連接度的統(tǒng)計值，k就是連接度值高蜂，每個power值提供了max, median, max3種連接度的統(tǒng)計量
# 對連接度的均值進行可視化
# 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")
WGCNA-4
★
★ 橫坐標為軟閾值的梯度备恤，第一幅圖的縱坐標為無標度網(wǎng)絡(luò)適應(yīng)系數(shù)稿饰，越大越好；第二幅圖的縱坐標為節(jié)點的平均連通度露泊，越小越好
2.1.2 一步構(gòu)建網(wǎng)絡(luò)與模塊檢測
# 打開多線程
enableWGCNAThreads(nThreads = 10)
# LUAD_Expr表達數(shù)據(jù)喉镰，TOMType拓撲重疊矩陣計算方式，minModuleSize用于模塊檢測的最小模塊尺寸,
# reassignThreshold 是否在模塊之間重新分配基因的p值比率閾值惭笑，mergeCutHeight 樹狀圖切割高度
# numericLabels 返回的模塊應(yīng)該用顏色（FALSE）還是數(shù)字（TRUE）標記,pamRespectsDendro樹狀圖相關(guān)參數(shù)
# saveTOMs 字符串的向量侣姆，saveTOMFileBase 包含包含共識拓撲重疊文件的文件名庫的字符串
net = blockwiseModules(LUAD_Expr, power = sft$powerEstimate,
maxBlockSize = 2000, # 最大模塊數(shù)量
TOMType = "unsigned", minModuleSize = 30, reassignThreshold = 0, # 構(gòu)建無尺度網(wǎng)絡(luò)，最小模塊基因數(shù)為30
mergeCutHeight = 0.25, # 需要合并模塊的閾值
numericLabels = TRUE, # 以數(shù)字作為模塊的名字
pamRespectsDendro = FALSE, saveTOMs = TRUE,
saveTOMFileBase = "LUADTOM", verbose = 3)
# 回到網(wǎng)絡(luò)分析沉噩，查看標識了多少個模塊以及模塊大小
table(net$colors)
##
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56
## 4348 9126 1444 1202 1158 1141 969 725 581 531 444 425 391 359 324 317 280 274 263 252 225 186 178 152 150 142 137 128 125 118 115 115 114 112 112 108 108 106 100 96 96 92 92 90 89 88 83 82 81 81 81 80 80 79 78 78 77
## 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104
## 74 74 74 68 67 64 64 64 62 61 60 59 58 57 56 55 55 54 52 52 52 52 50 50 50 49 49 49 47 45 43 43 39 39 38 38 38 36 36 36 34 34 33 33 33 31 31 30
# 將標簽轉(zhuǎn)化為繪圖顏色
moduleColors = labels2colors(net$colors)
# 繪制樹狀圖和下面的模塊顏色
# dendroLabels樹狀圖標簽捺宗。設(shè)置為FALSE完全禁用樹狀圖標簽；設(shè)置為NULL使用的行標簽LUAD_Expr
# addGuide是否應(yīng)在樹狀圖中添加垂直的“指導線”川蒙？線條使識別單個樣本的顏色代碼更加容易蚜厉。
plotDendroAndColors(net$dendrograms[[1]], moduleColors[net$blockGenes[[1]]], "Module colors",
dendroLabels = FALSE, hang = 0.03, addGuide = TRUE, guideHang = 0.05)
WGCNA-5
★ 指示有104個模塊，按大小降序標記為1至104畜眨，大小范圍為9126至30個基因弯囊； 標簽0保留用于所有模塊外部的基因
★ 無法聚類到模塊中的基因會標示為灰色痰哨，如果灰色區(qū)域較多，可能由于樣本中基因共表達趨勢不明顯匾嘱，可能需要調(diào)整基因過濾的方法
# 保存后續(xù)分析所需的模塊分配和模塊本征信息斤斧，可由recutBlockwiseTrees函數(shù)應(yīng)用修改后的條件而不必重新計算網(wǎng)絡(luò)和聚類樹狀圖
moduleLabels = net$colors
moduleColors = labels2colors(net$colors)
MEs = net$MEs;
geneTree = net$dendrograms[[1]];
save(MEs, moduleLabels, moduleColors, geneTree, file = "TCGA-LUAD-02-networkConstruction-auto.RData")
2.2 分步網(wǎng)絡(luò)構(gòu)建和模塊檢測
2.3 逐塊網(wǎng)絡(luò)構(gòu)建和模塊檢測
3. 分步網(wǎng)絡(luò)構(gòu)建和模塊檢測
3.1 量化模塊-特質(zhì)關(guān)聯(lián)
在此分析中，我們想確定與所測量的臨床特征顯著相關(guān)的模塊霎烙。 由于我們已經(jīng)為每個模塊建立了一個概要文件（特征基因）撬讽，因此我們只需將特征基因與外部特征相關(guān)聯(lián)，然后尋找最重要的關(guān)聯(lián)悬垃，實際上是計算模塊的ME值與表型的相關(guān)系數(shù)：
# Define numbers of genes and samples
# 獲得基因數(shù)和樣本數(shù)
nGenes = ncol(LUAD_Expr);
nSamples = nrow(LUAD_Expr);
# Recalculate MEs with color labels
# 用彩色標簽重新計算MEs
# 在給定的單個數(shù)據(jù)集中計算模塊的模塊本征基因
MEs0 = moduleEigengenes(LUAD_Expr, moduleColors)$eigengenes
# 對給定的（特征）向量進行重新排序游昼，以使相似的向量（通過相關(guān)性度量）彼此相鄰
MEs = orderMEs(MEs0)
# 計算module的ME值與表型的相關(guān)系數(shù)（pearson）
moduleTraitCor = cor(MEs, data_Traits, use = "p");
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples);
# 使用的是其他相關(guān)性方法，則可以使用 bicorAndPvalue 函數(shù)來計算顯著性
# modTraitCorP = bicorAndPvalue(MEs_col, design)
# modTraitCor = modTraitCorP$bicor
# modTraitP = modTraitCorP$p
names(MEs)
## [1] "MEdarkgrey" "MEyellowgreen" "MEivory" "MEnavajowhite1" "MEindianred3" "MElightcyan" "MEgreenyellow" "MEskyblue1" "MEbisque4" "MEblack" "MEbrown" "MEmagenta" "MEmidnightblue" "MEpalevioletred3"
## [15] "MEskyblue4" "MEmediumpurple1" "MEthistle1" "MEblue2" "MEorangered3" "MEdarkviolet" "MEplum1" "MEthistle" "MEdarkred" "MEwhite" "MEfirebrick3" "MEblueviolet" "MEsienna4" "MElightsteelblue"
## [29] "MElightyellow" "MEpalevioletred2" "MEcoral3" "MEpink4" "MEhoneydew" "MEyellow3" "MEnavajowhite2" "MEcoral1" "MEtan" "MEblue" "MEdarkolivegreen" "MEthistle3" "MEdarkseagreen4" "MEplum3"
## [43] "MEbrown2" "MElavenderblush3" "MEhoneydew1" "MElightslateblue" "MEdarkolivegreen4" "MEdarkturquoise" "MEviolet" "MEcoral" "MEpink" "MEmagenta4" "MEfirebrick4" "MEpurple" "MEred" "MElightcoral"
## [57] "MEyellow" "MEmediumpurple4" "MEskyblue" "MEcyan" "MEdarkgreen" "MElightpink4" "MEskyblue2" "MEmediumorchid" "MElightsteelblue1" "MEmaroon" "MEmediumpurple2" "MEplum2" "MEsteelblue" "MEdarkorange2"
## [71] "MEplum" "MEdarkmagenta" "MEorange" "MEdarkorange" "MEsalmon" "MElightgreen" "MEturquoise" "MEroyalblue" "MEyellow4" "MEorangered4" "MEpaleturquoise" "MElightblue4" "MEindianred4" "MEbrown4"
## [85] "MEsalmon4" "MEsienna3" "MEthistle2" "MEgrey60" "MEgreen" "MEorangered1" "MEmediumpurple3" "MEdarkslateblue" "MElightpink3" "MEantiquewhite2" "MElightcyan1" "MEfloralwhite" "MElavenderblush2" "MEdarkolivegreen2"
## [99] "MEskyblue3" "MEsalmon2" "MEcoral2" "MEsaddlebrown" "MEantiquewhite4" "MEdarkseagreen3" "MEgrey"
# sizeGrWindow(10,6)
# 顯示相關(guān)性及其p值
textMatrix = paste(signif(moduleTraitCor, 2), "\n(",signif(moduleTraitPvalue, 1), ")", sep = "");
dim(textMatrix) = dim(moduleTraitCor)
# Display the correlation values within a heatmap plot\
# ySymbols 當ylabels使用時所使用的其他標簽尝蠕； colorLabels 應(yīng)該使用顏色標簽嗎
# colors 顏色烘豌； textMatrix 單元格名字
labeledHeatmap(Matrix = moduleTraitCor, xLabels = names(data_Traits), yLabels = names(MEs), ySymbols = names(MEs),
colorLabels = FALSE,colors = greenWhiteRed(50), textMatrix = textMatrix,setStdMargins = FALSE,
cex.text = 0.2, zlim = c(-1,1),
main = paste("Module-trait relationships"))
WGCNA-6
★ 計算出模塊與表型之間相關(guān)性之后，可以挑選最相關(guān)的那些模塊來進行后續(xù)分析看彼。但是廊佩，模塊本身可能還包含很多的基因，還需要進一步識別關(guān)鍵基因
3.2 基因與性狀和重要模塊的關(guān)系：基因重要性和模塊成員
量化陣列上所有基因與每個模塊的相似性尋找重要模塊:
# 定義包含數(shù)據(jù)特征權(quán)重列的變量權(quán)重
days_to_death = as.data.frame(data_Traits$days_to_death);
names(days_to_death) = "days_to_death"
geneModuleMembership = as.data.frame(cor(LUAD_Expr, MEs, use = "p"));
# 模塊的名稱(顏色)substring提取文本從第3個字母開始（此處輸入法有問題靖榕，不用"#")
# modNames = substring(names(MEs), 3)
# 基因和模塊的相關(guān)系數(shù)
geneModuleMembership = as.data.frame(cor(LUAD_Expr, 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="");
# gene和性狀的關(guān)系
geneTraitSignificance = as.data.frame(cor(LUAD_Expr, days_to_death, use = "p"));
GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples));
names(geneTraitSignificance) = paste("GS.", names(days_to_death), sep="");
names(GSPvalue) = paste("p.GS.", names(days_to_death), sep = "");
3.3 模內(nèi)分析：鑒定具有高GS和MM的基因
使用GS和MM度量标锄，我們可以鑒定出對“days_to_death”以及在感興趣的模塊中具有較高模塊成員性具有重要意義的基因。 例如茁计，我們看一下與“days_to_death”關(guān)聯(lián)“brown”模塊料皇。 我們在“brown”模塊中繪制了基因重要性與模塊成員關(guān)系的散點圖。 在此模塊中星压，GS和MM之間存在高度顯著的負相關(guān)性
# 模型顏色
module = "brown"
# 匹配列
column = match(module, modNames);
moduleGenes = moduleColors==module;
#sizeGrWindow(7, 7);
par(mfrow = c(1,1));
# MM <- abs(geneModuleMembership[moduleGenes, column])
# GS <- abs(geneTraitSignificance[moduleGenes, 1])
# 畫散點圖
verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
abs(geneTraitSignificance[moduleGenes, 1]),
xlab = paste("Module Membership in", module, "module"),
ylab = "Gene significance for body weight",
main = paste("Module membership vs. gene significance\n"),
cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
WGCNA-7
★ GS是Gene Signicance践剂，描述的是某一個基因與性狀的相關(guān)性
★ MM是Module Membership，描述的是某一個基因與模塊之間的相關(guān)性
3.4 導出模塊網(wǎng)絡(luò)
file <- "~/TCGA-LUAD.net"
TOM <- TOMsimilarityFromExpr(LUAD_Expr, power = 2, networkType = "unsigned")
dimnames(TOM) <- list(colnames(LUAD_Expr), colnames(LUAD_Expr))
modTOM <- TOM[moduleGenes, moduleGenes]
cyt <- exportNetworkToCytoscape(
modTOM,
edgeFile = paste0(file, module, ".edges.txt"),
nodeFile = paste0(file, module, ".nodes.txt"),
weighted = TRUE,
threshold = 0, # threshold 默認為0.5, 可以根據(jù)自己的需要調(diào)整
nodeNames = moduleGenes,
nodeAttr = module
)
4. 使用WGCNA進行網(wǎng)絡(luò)可視化
4.1 顯示基因網(wǎng)絡(luò)
可視化加權(quán)網(wǎng)絡(luò)的一種方法是繪制其熱圖娜膘，熱圖的每一行和每一列都對應(yīng)一個基因舷手。 熱圖可以描述鄰接或拓撲重疊，淺色表示低鄰接（重疊）劲绪，而深色表示更高的鄰接（重疊）。 另外盆赤，沿著熱圖的頂部和左側(cè)繪制了基因樹狀圖和模塊顏色贾富。
# 重新計算拓撲重疊：通過保存TOM可以更有效地完成此操作
dissTOM = 1-TOMsimilarityFromExpr(LUAD_Expr, power = 2);
# 變換dissTOM（對dissTOM矩陣進行指數(shù)轉(zhuǎn)換，使中等強度的關(guān)系更容易在熱圖上展示出來）牺六；
plotTOM = dissTOM^7;
# #將對角線數(shù)據(jù)設(shè)為NA
diag(plotTOM) = NA;
# Call the plot function
# sizeGrWindow(9,9)
# 基因的聚類樹聚類時的距離為1-TOM值結(jié)合基因間的距離颤枪，即1-TOM值，用熱圖展示
# TOMplot(plotTOM, geneTree, moduleColors, main = "Network heatmap plot, all genes") # 此處我的基因數(shù)量過多淑际？會報錯：Error in x[, iy] : subscript out of bounds
★ 數(shù)據(jù)不合理畏纲，可根據(jù)探針集或基因的平均表達量或方差（如中位數(shù)或絕對中位差）重新進行過濾扇住，低表達或無變化的基因通常代表噪音。
# 限制基因數(shù)量以加快繪圖速度盗胀。 注意艘蹋，一個基因子集的基因樹狀圖通常看起來與所有基因的基因樹狀圖不同：
nSelect = 400
# 為了反復(fù)重現(xiàn)結(jié)果, 這里設(shè)置了隨機種子票灰；
set.seed(10);
select = sample(nGenes, size = nSelect);
# 提取抽取基因相應(yīng)的TOM矩陣
selectTOM = dissTOM[select, select];
# 重新畫聚類圖
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;
# 更換配色
# mycolor <- gplots::colorpanel(250, 'red', "orange", 'lemonchiffon')
TOMplot(plotDiss, selectTree, selectColors,
main = "Network heatmap plot, selected genes",
# col = mycolor)
WGCNA-8
4.2 可視化特征基因網(wǎng)絡(luò)
研究找到的模塊之間的關(guān)系通常很有趣女阀。 可以使用特征基因作為代表特征，并通過特征基因相關(guān)性來量化模塊相似性屑迂。 該軟件包包含一個方便的函數(shù)plotEigengeneNetworks浸策，該函數(shù)生成特征基因網(wǎng)絡(luò)的摘要圖。 通常惹盼，向特征基因添加臨床特征（或多個特征）以了解特征如何適合特征基因網(wǎng)絡(luò)是有益的
# Recalculate module eigengenes
# 重新計算基因特征值
MEs = moduleEigengenes(LUAD_Expr, moduleColors)$eigengenes
# Isolate weight from the clinical traits
days_to_death = as.data.frame(data_Traits$days_to_death);
names(days_to_death) = "days_to_death"
# Add the weight to existing module eigengenes
MET = orderMEs(cbind(MEs, days_to_death))
# Plot the relationships among the eigengenes and the trait
# sizeGrWindow(5,7.5);
par(cex = 0.2)
# 畫樹形圖
# marDendro給出樹狀圖的邊距設(shè)置庸汗，marHeatmap熱圖邊距設(shè)置
plotEigengeneNetworks(MET, "", marDendro = c(0,4,1,2), marHeatmap = c(3,4,1,2), cex.lab = 0.8, yLabelsAngle = 90)
WGCNA-9
拆分樹狀圖和熱圖圖
# Plot the dendrogram
# sizeGrWindow(6,6);
par(cex = 1.0)
plotEigengeneNetworks(MET, "Eigengene dendrogram", marDendro = c(0,4,2,0),
plotHeatmaps = FALSE)
# Plot the heatmap matrix (note: this plot will overwrite the dendrogram plot)
par(cex = 1.0)
plotEigengeneNetworks(MET, "Eigengene adjacency heatmap", marHeatmap = c(3,4,2,2),plotDendrograms = FALSE, yLabelsAngle = 90)
