在做單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析的時候，獲得分群后我們希望知道每個群的marker基因是什么以此來描述細(xì)胞群的特征。其技術(shù)本質(zhì)就是利用某種規(guī)則獲得樣本（高変基因）或者每個亞群（差異基因）的一個基因list。WGCNA為我們提供了一套有別于高変基因、差異分析的規(guī)則，找到某些細(xì)胞中有關(guān)聯(lián)作用的基因list（他們稱作模塊）。它與只挑選差異基因相比，WGCNA可以從成千上萬的基因中挑選出高度相關(guān)的基因的簇（模塊），并將模塊與外部樣本性狀關(guān)聯(lián)，找出與樣本性狀高度相關(guān)的模塊。然后就可以進行模塊內(nèi)分析。

在一般的描述中：

To systematically investigate the genetic program dynamics, we performed weighted gene co-expression network analysis (WGCNA) on 2,464 genes that were variably expressedin trophoblast cells between different developmental stages. WGCNA identified 8 gene modules, each of which contains a set of genes that tend to be co-expressed at a certain development stage.

在漢語世界WGCNA已經(jīng)有大量的教程了（大多是基于bulk-RNA或者是探針數(shù)據(jù)），這是一種相關(guān)性的技術(shù)工具，可以說把相關(guān)性應(yīng)用到了新的高度。加權(quán)基因共表達網(wǎng)絡(luò)分析 (WGCNA, Weighted correlation network
analysis)是用來描述不同樣品（我們的cell-barcode）之間基因關(guān)聯(lián)模式的系統(tǒng)生物學(xué)方法，可以用來鑒定高度協(xié)同變化的基因集，并根據(jù)基因集的內(nèi)連性和基因集與表型之間的關(guān)聯(lián)鑒定marker gene 或治療靶點。

WGCNA 分析一般需要一個表達矩陣，和一個表型矩陣（如果有的話），所以作為一個工具箱，是可以應(yīng)用在單細(xì)胞中的，但是鑒于單細(xì)胞不同技術(shù)的數(shù)據(jù)特點不同，在用的時候我們需要注意一下。這不是一個新鮮的問題，在傳統(tǒng)的WGCNA中就有這個問題：關(guān)鍵問題答疑：WGCNA的輸入矩陣到底是什么格式？

大佬的推薦：單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析的時候可以加上wgcna

什么是WGCNA

首先我們還是要知道一下WGCNA的基本概念。在我們想要學(xué)習(xí)WGCNA的時候，或者干脆很多人就是從一文看懂WGCNA 分析(2019更新版)中了解和學(xué)習(xí)WGCNA的。關(guān)于這分析我們需要明白兩個核心的概念：相關(guān)性、聚類。這兩個都不是陌生的概念，相關(guān)性就是兩個變量之間的協(xié)同變化，同增同減還是向相而行？聚類就是計算兩個變量之間的距離，物以類聚人以群分。

要達到這種技術(shù)目標(biāo)，需要一套新的概念：

• Co-expression network(共表達網(wǎng)絡(luò))：共表達網(wǎng)絡(luò)定義為無向的、加權(quán)的基因網(wǎng)絡(luò)。這樣一個網(wǎng)絡(luò)的節(jié)點對應(yīng)于基因，基因之間的邊代表基因表達量的相關(guān)性，加權(quán)是將相關(guān)性的絕對值提高到冪β≥1(軟閾值)，加權(quán)基因共表達網(wǎng)絡(luò)的構(gòu)建以犧牲低相關(guān)性為代價，強調(diào)高相關(guān)性。
• Module(模塊) ： 模塊是高度互連的基因簇。模塊對應(yīng)于正相關(guān)的基因。這里的加權(quán)的網(wǎng)絡(luò)就等于鄰接矩陣。通過冪鄰接轉(zhuǎn)換，就強化了高相關(guān)性基因的關(guān)系，弱化了相關(guān)性基因的關(guān)系。
• Connectivity(連接度)： 對于每個基因，連接性(也稱為度)被定義為與其他基因的連接強度之和：在共表達網(wǎng)絡(luò)中，連接度衡量一個基因與所有其他網(wǎng)絡(luò)基因的相關(guān)性。
• Intramodular connectivity(模塊內(nèi)連接度)： 模塊內(nèi)鏈接度衡量給定基因相對于特定模塊的基因的連接或共表達程度。模塊內(nèi)連接度可以做為Module membership的度量。
• Module eigengene E： 給定模塊的第一主成分，代表整個模塊的基因表達譜
• Module Membership(MM)： 對于每個基因，我們通過將其基因表達譜與模塊的Module eigengene相關(guān)性來定義Module Membership。
• hub gene ： 高度連接基因的縮寫，根據(jù)定義，它是共表達網(wǎng)絡(luò)模塊內(nèi)具有高連接度的基因。
• Gene significance(GS) ：
• 模塊的顯著性（module significance,Ms）： 定義為模塊包含的所有基因顯著性性的平均值，然后比較MS，一般MS越高，說明這個模塊與疾病之間的關(guān)聯(lián)度越高 。

下面我們來看一下WGCNA的官方流程：

• 構(gòu)建基因共表達網(wǎng)絡(luò)：使用加權(quán)的表達相關(guān)性。
• 識別基因集：基于加權(quán)相關(guān)性，進行層級聚類分析，并根據(jù)設(shè)定標(biāo)準(zhǔn)切分聚類結(jié)果，獲得不同的基因模塊，用聚類樹的分枝和不同顏色表示。
• 如果有表型信息，計算基因模塊與表型的相關(guān)性，鑒定性狀相關(guān)的模塊。
• 研究模塊之間的關(guān)系，從系統(tǒng)層面查看不同模塊的互作網(wǎng)絡(luò)。
• 從關(guān)鍵模塊中選擇感興趣的驅(qū)動基因，或根據(jù)模塊中已知基因的功能推測未知基因的功能。
• 導(dǎo)出TOM矩陣，繪制相關(guān)性網(wǎng)絡(luò)圖。

看了這個大致的流程我們有了一種感覺：這像一種特征選擇器，根據(jù)相關(guān)性選基因模塊。那么，當(dāng)我做WGCNA時，我在做什么?

• 獲得樣本（或亞群）的基因模塊
• 模塊與表型之間的關(guān)系
• 選出來的模塊的差異也是樣本（或亞群）的異質(zhì)性的表現(xiàn)

sc-WGCNA

在github上面還真有人創(chuàng)建了一個項目：https://github.com/milescsmith/scWGCNA

當(dāng)我們在單細(xì)胞數(shù)據(jù)分析中應(yīng)用WGCNA的時候第一個問題就還是：到底輸入矩陣是什么？ 我們知道在10X的scrna分析中Seurat有三個數(shù)據(jù)：

• count ： 原始count
• data ： 均一化之后
• scale.data： 標(biāo)準(zhǔn)化之后

在我們閱讀完WGCNA的一些教程之后，我們發(fā)現(xiàn)他需要的其實是均一化之后的數(shù)據(jù)，也就是一般的在Seurat對象的pbmc_small@assays$RNA@data·中的數(shù)據(jù)。表型數(shù)據(jù)一般就是pbmc_small@meta.data`。

在WGCNA的文檔中我們也常常看到這樣的建議：

不建議對少于15個樣本的數(shù)據(jù)集嘗試WGCNA。與其他分析方法一樣，更多的樣品通常會導(dǎo)致更可靠和更精確的結(jié)果。

這一點高通量的單細(xì)胞轉(zhuǎn)錄組一般是不怕的，細(xì)胞數(shù)一般在幾千。但是這也是一個挑戰(zhàn)：A）數(shù)據(jù)緯度高了計算量大；B)緯度高數(shù)據(jù)稀疏，相關(guān)性差，找不到明顯的模塊。

于是，我們看到Pseudocell的概念：

Tosches, M. A. et al. Evolution of pallium, hippocampus, and cortical cell types revealed by single-cell transcriptomics in reptiles. Science 360, 881-888

在用單細(xì)胞數(shù)據(jù)的WGCNA分析之前也是每個cluster隨機選一部分細(xì)胞構(gòu)成Pseudocell（局部bulk的方法）。怪不得我用原始的count矩陣做WGCNA的結(jié)果這么差呢。

傳統(tǒng)地，探針集或基因可以通過均值、絕對中位差(MAD)或方差進行過濾，因為低表達或不變的基因通常代表噪聲。用均值表達還是方差過濾是否更好尚有爭議，兩者都有優(yōu)缺點。不建議通過差異分析過濾掉基因。我們知道單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)的一個特點就是緯度高，數(shù)據(jù)稀疏，除了要考慮細(xì)胞的特殊處理之外，我們還可以過濾基因，如只用高変基因（FindVariableFeatures()）。這里請注意，也是不推薦只選用某一群的差異基因做的，因為某一群的差異基因，已經(jīng)是一個明顯的模塊了，這樣做很可能只得到很少的模塊。

在解決了數(shù)據(jù)數(shù)據(jù)的處理之后，我們就可以用豐富的WGCNA教程來分析我們的單細(xì)胞數(shù)據(jù)了。

WGCNA in action

在討論了一般的規(guī)則之后，我們就可以著手跑自己的單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)集了。再次明確一下數(shù)據(jù)集的一些建議：

• 基因過濾，可以挑選一部分基因做
• 細(xì)胞過濾，A選擇某一群細(xì)胞看某一群的基因表達模塊；B整個樣本做；C如果細(xì)胞數(shù)據(jù)比較離散可以考慮我們上面提到的構(gòu)造Pseudocell再做
• 數(shù)據(jù)一般使用均一化之后的

下面我們就參考一文看懂WGCNA 分析(2019更新版)來做一遍。

載入R包讀取數(shù)據(jù)：

library(WGCNA)
library(Seurat)
library(tidyverse)
library(reshape2)
library(stringr)

pbmc <- readRDS('G:\\Desktop\\Desktop\\RStudio\\single_cell\\filtered_gene_bc_matrices\\hg19pbmc_tutorial.rds')

pbmc

An object of class Seurat 
13714 features across 2638 samples within 1 assay 
Active assay: RNA (13714 features)
 3 dimensional reductions calculated: pca, umap, tsne

 head(pbmc@meta.data)  #表型數(shù)據(jù)也有了
               orig.ident nCount_RNA nFeature_RNA percent.mt RNA_snn_res.0.5 seurat_clusters
AAACATACAACCAC     pbmc3k       2419          779  3.0177759               1               1
AAACATTGAGCTAC     pbmc3k       4903         1352  3.7935958               3               3
AAACATTGATCAGC     pbmc3k       3147         1129  0.8897363               1               1
AAACCGTGCTTCCG     pbmc3k       2639          960  1.7430845               2               2
AAACCGTGTATGCG     pbmc3k        980          521  1.2244898               6               6
AAACGCACTGGTAC     pbmc3k       2163          781  1.6643551               1               1

大家看到了，我們用的是10X的數(shù)據(jù)，細(xì)胞數(shù)很多，遠大于15，矩陣也以稀疏著稱，所以我們還是先把細(xì)胞捏一下。

datadf <- as.matrix(pbmc@assays$RNA@data )
idd1 <- pbmc@meta.data
Inter.id1<-cbind(rownames(idd1),idd1$seurat_clusters)
rownames(Inter.id1)<-rownames(idd1)
colnames(Inter.id1)<-c("CellID","Celltype")
Inter.id1<-as.data.frame(Inter.id1)
head(Inter.id1)
Inter1<-datadf[,Inter.id1$CellID]
Inter2<-as.matrix(Inter1)
Inter2[1:4,1:4]

pseudocell.size = 10 ## 10 test
new_ids_list1 = list()
length(levels(Inter.id1$Celltype))

for (i in 1:length(levels(Inter.id1$Celltype))) {
  cluster_id = levels(Inter.id1$Celltype)[i]
  cluster_cells <- rownames(Inter.id1[Inter.id1$Celltype == cluster_id,])
  cluster_size <- length(cluster_cells)     
  pseudo_ids <- floor(seq_along(cluster_cells)/pseudocell.size)
  pseudo_ids <- paste0(cluster_id, "_Cell", pseudo_ids)
  names(pseudo_ids) <- sample(cluster_cells)    
  new_ids_list1[[i]] <- pseudo_ids      
}

new_ids <- unlist(new_ids_list1)
new_ids <- as.data.frame(new_ids)
head(new_ids)
new_ids_length <- table(new_ids)
new_ids_length

new_colnames <- rownames(new_ids)  ###add
#rm(all.data1)
gc()
colnames(datadf)  
all.data<-datadf[,as.character(new_colnames)] ###add
all.data <- t(all.data)###add
new.data<-aggregate(list(all.data[,1:length(all.data[1,])]),
                    list(name=new_ids[,1]),FUN=mean)
rownames(new.data)<-new.data$name
new.data<-new.data[,-1]
new_ids_length<-as.matrix(new_ids_length)##

short<-which(new_ids_length<10)##
new_good_ids<-as.matrix(new_ids_length[-short,])##
result<-t(new.data)[,rownames(new_good_ids)]
dim(result)

13714   252

還剩下252個細(xì)胞，我們再把基因過濾一下。

pbmc <- FindVariableFeatures(pbmc,nfeatures = 5000)
colnames(result)[grepl("[12]_Cel",colnames(result))]
Cluster1 <- result[intersect(Seurat::VariableFeatures(pbmc),rownames(result)),]

最終我們的表達譜是這樣的：

 Cluster1[1:4,1:4]
         1_Cell1  1_Cell10  1_Cell11 1_Cell12
PPBP   0.0000000 0.0000000 0.0000000 0.000000
LYZ    1.4139673 0.4601388 0.3870654 1.123586
S100A9 0.1905211 0.1616515 0.0000000 0.507018
IGLL5  0.0000000 0.0000000 0.0000000 0.000000

 dim(Cluster1)
[1] 4273  252

處理好表達譜，下面的流程就常見了。

WGCNA基本參數(shù)設(shè)置：

type = "unsigned"  # 官方推薦 "signed" 或 "signed hybrid"
corType = "pearson" # 相關(guān)性計算  官方推薦 biweight mid-correlation & bicor  corType: pearson or bicor 
corFnc = ifelse(corType=="pearson", cor, bicor)
corFnc
maxPOutliers = ifelse(corType=="pearson",1,0.05) # 對二元變量，如樣本性狀信息計算相關(guān)性時， # 或基因表達嚴(yán)重依賴于疾病狀態(tài)時，需設(shè)置下面參數(shù)
# 關(guān)聯(lián)樣品性狀的二元變量時，設(shè)置
robustY = ifelse(corType=="pearson",T,F)
dataExpr  <- as.matrix(Cluster1)

根據(jù)表達量再做一次篩選。

## 篩選中位絕對偏差前75%的基因，至少MAD大于0.01
## 篩選后會降低運算量，也會失去部分信息
## 也可不做篩選，使MAD大于0即可

m.mad <- apply(dataExpr,1,mad)
dataExprVar <- dataExpr[which(m.mad > 
                                max(quantile(m.mad, probs=seq(0, 1, 0.25))[2],0.01)),]

## 轉(zhuǎn)換為樣品在行，基因在列的矩陣
dataExpr <- as.data.frame(t(dataExprVar))
dim(dataExpr)
head(dataExpr)[,1:8]
               LYZ    S100A9      GNLY      FTL     FTH1    S100A8     CD74      NKG7
1_Cell1  1.4139673 0.1905211 0.1573319 2.786124 3.193984 0.0000000 2.456273 0.4739367
1_Cell10 0.4601388 0.1616515 0.2451389 2.506929 3.112310 0.1558935 1.446417 0.6014867
1_Cell11 0.3870654 0.0000000 0.3857747 2.448856 3.327483 0.0000000 1.460063 0.2179642
1_Cell12 1.1235861 0.5070180 0.0000000 2.965330 3.084265 0.0000000 2.084446 1.0519813
1_Cell13 1.0280862 0.1635541 0.0000000 2.969338 3.044814 0.0000000 1.554540 0.2370180
1_Cell14 1.0568580 0.1567888 0.0000000 2.962782 3.009191 0.0000000 1.658674 0.3956653

## 檢測缺失值
gsg = goodSamplesGenes(dataExpr, verbose = 3)
gsg$allOK
gsg$goodSamples

if (!gsg$allOK){
  # Optionally, print the gene and sample names that were removed:
  if (sum(!gsg$goodGenes)>0) 
    printFlush(paste("Removing genes:", 
                     paste(names(dataExpr)[!gsg$goodGenes], collapse = ",")));
  if (sum(!gsg$goodSamples)>0) 
    printFlush(paste("Removing samples:", 
                     paste(rownames(dataExpr)[!gsg$goodSamples], collapse = ",")));
  # Remove the offending genes and samples from the data:
  dataExpr = dataExpr[gsg$goodSamples, gsg$goodGenes]
}

nGenes = ncol(dataExpr)
nSamples = nrow(dataExpr)

dim(dataExpr)
[1]  252 1049

注意一下，在篩選的時候，有些基因是要保留的，哪些基因呢？就是那些你關(guān)注的基因。

## 查看是否有離群樣品
sampleTree = hclust(dist(dataExpr), method = "average")
plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="")

這要是幾千個細(xì)胞的話，簡直一團黑。

powers = c(c(1:10), seq(from = 12, to=30, by=2))
sft = pickSoftThreshold(dataExpr, powerVector=powers, 
                        networkType="signed", verbose=5)

關(guān)鍵就是理解pickSoftThreshold函數(shù)及其返回的對象，最佳的beta值就是sft$powerEstimate

par(mfrow = c(1,2))
cex1 = 0.9
# 橫軸是Soft threshold (power)，縱軸是無標(biāo)度網(wǎng)絡(luò)的評估參數(shù)，數(shù)值越高，
# 網(wǎng)絡(luò)越符合無標(biāo)度特征 (non-scale)
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     xlab="Soft Threshold (power)",
     ylab="Scale Free Topology Model Fit,signed R^2",type="n",
     main = paste("Scale independence"))
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     labels=powers,cex=cex1,col="red")
# 篩選標(biāo)準(zhǔn)。R-square=0.85
abline(h=0.85,col="red")

# Soft threshold與平均連通性
plot(sft$fitIndices[,1], sft$fitIndices[,5],
     xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
     main = paste("Mean connectivity"))
text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, 
     cex=cex1, col="red")

power = sft$powerEstimate
softPower  = power
softPower

5

參數(shù)beta取值默認(rèn)是1到30，上述圖形的橫軸均代表權(quán)重參數(shù)β，左圖縱軸代表對應(yīng)的網(wǎng)絡(luò)中l(wèi)og(k)與log(p(k))相關(guān)系數(shù)的平方。相關(guān)系數(shù)的平方越高，說明該網(wǎng)絡(luò)越逼近無網(wǎng)路尺度的分布。右圖的縱軸代表對應(yīng)的基因模塊中所有基因鄰接函數(shù)的均值。最佳的beta值就是sft$powerEstimate，已經(jīng)被保存到變量了，不需要知道具體是什么，后面的代碼都用這個即可，在本例子里面是5。

即使你不理解它，也可以使用代碼拿到合適“軟閥值（soft thresholding power）”beta進行后續(xù)分析。

# 無向網(wǎng)絡(luò)在power小于15或有向網(wǎng)絡(luò)power小于30內(nèi)，沒有一個power值可以使
# 無標(biāo)度網(wǎng)絡(luò)圖譜結(jié)構(gòu)R^2達到0.8，平均連接度較高如在100以上，可能是由于
# 部分樣品與其他樣品差別太大。這可能由批次效應(yīng)、樣品異質(zhì)性或?qū)嶒灄l件對
# 表達影響太大等造成?？梢酝ㄟ^繪制樣品聚類查看分組信息和有無異常樣品。
# 如果這確實是由有意義的生物變化引起的，也可以使用下面的經(jīng)驗power值。
if (is.na(power)){
  power = ifelse(nSamples<20, ifelse(type == "unsigned", 9, 18),
                 ifelse(nSamples<30, ifelse(type == "unsigned", 8, 16),
                        ifelse(nSamples<40, ifelse(type == "unsigned", 7, 14),
                               ifelse(type == "unsigned", 6, 12))       
                 )
  )
}

power

關(guān)鍵一步：

#一步法網(wǎng)絡(luò)構(gòu)建：One-step network construction and module detection##
# power: 上一步計算的軟閾值
# maxBlockSize: 計算機能處理的最大模塊的基因數(shù)量 (默認(rèn)5000)；
#  4G內(nèi)存電腦可處理8000-10000個，16G內(nèi)存電腦可以處理2萬個，32G內(nèi)存電腦可
#  以處理3萬個
#  計算資源允許的情況下最好放在一個block里面。
# corType: pearson or bicor
# numericLabels: 返回數(shù)字而不是顏色作為模塊的名字，后面可以再轉(zhuǎn)換為顏色
# saveTOMs：最耗費時間的計算，存儲起來，供后續(xù)使用
# mergeCutHeight: 合并模塊的閾值，越大模塊越少

# type = unsigned

cor <- WGCNA::cor

net = blockwiseModules(dataExpr, power = power, maxBlockSize = nGenes,#nGenes
                       TOMType = "unsigned", minModuleSize = 10,
                       reassignThreshold = 0, mergeCutHeight = 0.25,
                       numericLabels = TRUE, pamRespectsDendro = FALSE,
                       saveTOMs=TRUE, corType = corType, 
                       maxPOutliers=maxPOutliers, loadTOMs=TRUE,
                       saveTOMFileBase = paste0("dataExpr", ".tom"),
                       verbose = 3)

計算過程：

 Calculating module eigengenes block-wise from all genes
   Flagging genes and samples with too many missing values...
    ..step 1
 ..Working on block 1 .
    TOM calculation: adjacency..
    ..will not use multithreading.
     Fraction of slow calculations: 0.000000
    ..connectivity..
    ..matrix multiplication (system BLAS)..
    ..normalization..
    ..done.
   ..saving TOM for block 1 into file dataExpr.tom-block.1.RData
 ....clustering..
 ....detecting modules..
 ....calculating module eigengenes..
 ....checking kME in modules..
     ..removing 109 genes from module 1 because their KME is too low.
     ..removing 17 genes from module 2 because their KME is too low.
 ..merging modules that are too close..
     mergeCloseModules: Merging modules whose distance is less than 0.25
       Calculating new MEs...

我們肯定對這個返回的對象比較感興趣啦,，挨個查看一下。

table(net$colors)
net$unmergedColors
head(net$MEs)
net$goodSamples
net$goodGenes
net$TOMFiles
net$blockGenes
net$blocks
net$MEsOK
head(net$MEs)
table(net$colors)

 0   1   2   3   4   5 
546 291 111  55  29  17 

這里用不同的顏色來代表那些所有的模塊，其中灰色默認(rèn)是無法歸類于任何模塊的那些基因，如果灰色模塊里面的基因太多，那么前期對表達矩陣挑選基因的步驟可能就不太合適。

## 灰色的為**未分類**到模塊的基因。
# Convert labels to colors for plotting
moduleLabels = net$colors
moduleColors = labels2colors(moduleLabels)
moduleColors
# Plot the dendrogram and the module colors underneath
# 如果對結(jié)果不滿意，還可以recutBlockwiseTrees，節(jié)省計算時間
plotDendroAndColors(net$dendrograms[[1]], moduleColors[net$blockGenes[[1]]],
                    "Module colors",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)

我們這個至少不是全灰的啊。

# module eigengene, 可以繪制線圖，作為每個模塊的基因表達趨勢的展示
MEs = net$MEs

### 不需要重新計算，改下列名字就好
### 官方教程是重新計算的，起始可以不用這么麻煩
MEs_col = MEs
colnames(MEs_col) = paste0("ME", labels2colors(
  as.numeric(str_replace_all(colnames(MEs),"ME",""))))
MEs_col = orderMEs(MEs_col)

# 根據(jù)基因間表達量進行聚類所得到的各模塊間的相關(guān)性圖
# marDendro/marHeatmap 設(shè)置下、左、上、右的邊距
head(MEs_col)
?plotEigengeneNetworks
plotEigengeneNetworks(MEs, "Eigengene adjacency heatmap", 
                      marDendro = c(3,3,2,4),
                      marHeatmap = c(3,4,2,2),
                      plotDendrograms = T,
                      xLabelsAngle = 90)

計算每個模塊的特征向量基因，為某一特定模塊第一主成分基因E。

 計算每個模塊的特征向量基因，為某一特定模塊第一主成分基因E。代表了該模塊內(nèi)基因表達的整體水平
MEList = moduleEigengenes(dataExpr, colors = dynamicColors)
MEs = MEList$eigengenes
# 計算根據(jù)模塊特征向量基因計算模塊相異度：
MEDiss = 1-cor(MEs);
# Cluster module eigengenes
METree = hclust(as.dist(MEDiss), method = "average");
# Plot the result

plotEigengeneNetworks(MEs, 
                      "Eigengene adjacency heatmap", 
                      marHeatmap = c(3,4,2,2), 
                      plotDendrograms = FALSE, 
                      xLabelsAngle = 90) 

畫出指定模塊表達量的熱圖：

which.module="turquoise"; 
ME=mergedMEs[, paste("ME",which.module, sep="")]
par(mfrow=c(2,1), mar=c(0,4.1,4,2.05))
plotMat(t(scale(dataExpr[,moduleColors==which.module ]) ),
        nrgcols=30,rlabels=F,rcols=which.module,
        main=which.module, cex.main=2)
par(mar=c(2,2.3,0.5,0.8))
barplot(ME, col=which.module, main="", cex.main=2,
        ylab="eigengene expression",xlab="array sample")

所有基因模塊關(guān)系：

load(net$TOMFiles, verbose=T)

## Loading objects:
##   TOM

TOM <- as.matrix(TOM)
TOM[1:4,1:4]
#dim(TOM2)

dissTOM = 1-TOM
# 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
table(moduleColors)
# 這一部分特別耗時，行列同時做層級聚類
TOMplot(plotTOM, net$dendrograms[[1]], moduleColors[net$blockGenes[[1]]], 
        main = "Network heatmap plot, all genes")

由于細(xì)胞名稱被我們改變了，原來的對象我們要構(gòu)造表型數(shù)據(jù)：

Cluster1[1:4,1:4]
mypbmc<- CreateSeuratObject(Cluster1)

mypbmc
mypbmc[["percent.mt"]] <- PercentageFeatureSet(mypbmc, pattern = "^CD")

mypbmc <- FindVariableFeatures(mypbmc, selection.method = "vst", nfeatures = 2000)
mypbmc %>% NormalizeData( normalization.method = "LogNormalize", scale.factor = 10000)%>%
  FindVariableFeatures( selection.method = "vst", nfeatures = 2000) %>%
  ScaleData(features=VariableFeatures(mypbmc),vars.to.regress = "percent.mt")  %>%
  RunPCA(features = VariableFeatures(object = mypbmc)) %>%
  FindNeighbors( dims = 1:10) %>%
  FindClusters( resolution = 0.5) %>%
  BuildClusterTree() %>%
  RunUMAP( dims = 1:10)  -> mypbmc

head(mypbmc@meta.data)

head(mypbmc@meta.data)
         orig.ident nCount_RNA nFeature_RNA percent.mt RNA_snn_res.0.5 seurat_clusters
1_Cell1           1   477.9562         1030   3.535886               0               0
1_Cell10          1   508.4991         1094   3.107892               0               0
1_Cell11          1   487.7897         1067   3.413325               0               0
1_Cell12          1   464.1321          954   3.189983               0               0
1_Cell13          1   475.1926         1014   3.752087               0               0
1_Cell14          1   502.9780         1072   3.327189               0               0

通過模塊與各種表型的相關(guān)系數(shù)，可以很清楚的挑選自己感興趣的模塊進行下游分析了。

moduleTraitCor_noFP <- cor(mergedMEs, mypbmc@meta.data, use = "p");
moduleTraitPvalue_noFP = corPvalueStudent(moduleTraitCor_noFP, nSamples); 
textMatrix_noFP <- paste(signif(moduleTraitCor_noFP, 2), "\n(", signif(moduleTraitPvalue_noFP, 1), ")", sep = ""); 
par(mar = c(10, 8.5, 3, 3)); 
labeledHeatmap(Matrix = moduleTraitCor_noFP, 
               xLabels = names(mypbmc@meta.data), 
               yLabels = names(mergedMEs), 
               ySymbols = names(mergedMEs), 
               colorLabels = FALSE, 
               colors = blueWhiteRed(50), 
               textMatrix = textMatrix_noFP,
               setStdMargins = FALSE, 
               cex.text = 0.65, 
               zlim = c(-1,1), 
               main = paste("Module-trait relationships")) 

根據(jù)性狀與模塊特征向量基因的相關(guān)性及pvalue來挖掘與性狀相關(guān)的模塊

library(pheatmap)
cor_ADR <- signif(WGCNA::cor(mypbmc@meta.data,mergedMEs,use="p",method="pearson"),5)

p.values <- corPvalueStudent(cor_ADR,nSamples=nrow(mypbmc@meta.data))

pheatmap(cor_ADR,display_numbers = matrix(ifelse(p.values <= 0.01, "**", ifelse(p.values<= 0.05 ,"*"," ")), nrow(p.values)),fontsize=18)

根據(jù)基因網(wǎng)絡(luò)顯著性，也就是性狀與每個基因表達量相關(guān)性在各個模塊的均值作為該性狀在該模塊的顯著性，顯著性最大的那個模塊與該性狀最相關(guān)：

GS1 <- as.numeric(WGCNA::cor(mypbmc@meta.data[,3],dataExpr,use="p",method="pearson"))
# 顯著性是絕對值：
GeneSignificance <- abs(GS1)
length(GeneSignificance)
length(mergedColors)
mypbmc
dim(dataExpr)
dim(Cluster1)
dim(mypbmc@meta.data)

# 獲得該性狀在每個模塊中的顯著性：
ModuleSignificance <- tapply(GeneSignificance,mergedColors,mean,na.rm=T)
ModuleSignificance

     blue     brown     green      grey turquoise    yellow 
0.3817643 0.1335957 0.1147941 0.1291072 0.3386081 0.1243168 

尋找與該性狀相關(guān)的樞紐基因（hub genes）,首先計算基因的內(nèi)部連接度和模塊身份，內(nèi)部連接度衡量的是基因在模塊內(nèi)部的地位，而模塊身份表明基因?qū)儆谀膫€模塊。

# 計算每個基因模塊內(nèi)部連接度，也就是基因直接兩兩加權(quán)相關(guān)性。
ADJ1=abs(cor(dataExpr,use="p"))^softPower 
# 根據(jù)上面結(jié)果和基因所屬模塊信息獲得連接度：
# 整體連接度 kTotal，模塊內(nèi)部連接度：kWithin，kOut=kTotal-kWithin， kDiff=kIn-kOut=2*kIN-kTotal 
Alldegrees1=intramodularConnectivity(ADJ1, moduleColors) 
head(Alldegrees1)

          kTotal   kWithin     kOut     kDiff
LYZ    37.799195 35.682342 2.116853 33.565490
S100A9 35.938946 33.908114 2.030833 31.877281
GNLY    5.284786  4.399469 0.885317  3.514152
FTL    44.989900 41.637720 3.352181 38.285539
FTH1   43.791667 41.146166 2.645501 38.500665
S100A8 27.131010 25.587161 1.543848 24.043313

# 注意模塊內(nèi)基于特征向量基因連接度評估模塊內(nèi)其他基因： de ne a module eigengene-based connectivity measure for each gene as the correlation between a the gene expression and the module eigengene
# 如 brown 模塊內(nèi)：kM Ebrown(i) = cor(xi, MEbrown) ， xi is the gene expression pro le of gene i and M Ebrown is the module eigengene of the brown module
# 而 module membership 與內(nèi)部連接度不同。MM 衡量了基因在全局網(wǎng)絡(luò)中的位置。
datKME=signedKME(dataExpr, MEs, outputColumnName="MM.")
datKME[1:4,1:4]

# 注意模塊內(nèi)基于特征向量基因連接度評估模塊內(nèi)其他基因： de ne a module eigengene-based connectivity measure for each gene as the correlation between a the gene expression and the module eigengene
# 如 brown 模塊內(nèi)：kM Ebrown(i) = cor(xi, MEbrown) ， xi is the gene expression pro le of gene i and M Ebrown is the module eigengene of the brown module
# 而 module membership 與內(nèi)部連接度不同。MM 衡量了基因在全局網(wǎng)絡(luò)中的位置。
datKME=signedKME(dataExpr, MEs, outputColumnName="MM.")
datKME[1:4,1:4]

選擇特定模塊的基因：

table(moduleColors)
module = "yellow";
# Select module probes
probes = colnames(dataExpr) ## 我們例子里面的probe就是基因名
inModule = (moduleColors==module);
modProbes = probes[inModule]; 
modProbes

 [1] "GIMAP5"  "PPA1"    "CD2"     "AQP3"    "MAL"     "UBAC2"  
 [7] "WTAP"    "SNX3"    "SLFN5"   "ITM2A"   "SUCLG2"  "RWDD1"  
[13] "IL32"    "GIMAP7"  "TRADD"   "CLPP"    "IL7R"    "TOB1"   
[19] "GOLGA7"  "CISH"    "RARRES3" "NOSIP"   "TCF7"    "FXYD5"  
[25] "CUTA"    "CD3D"    "LDHB"    "CD3E"    "RAN"    

導(dǎo)出作圖文件，主要模塊里面的基因直接的相互作用關(guān)系信息可以導(dǎo)出到cytoscape,VisANT等網(wǎng)絡(luò)可視化軟件。

## 也是提取指定模塊的基因名
# Select the corresponding Topological Overlap
modTOM = TOM[inModule, inModule];
dimnames(modTOM) = list(modProbes, modProbes)
vis = exportNetworkToVisANT(modTOM,
                            file = paste("VisANTInput-", module, ".txt", sep=""),
                            weighted = TRUE,
                            threshold = 0)
?exportNetworkToCytoscape
cyt = exportNetworkToCytoscape(
  modTOM,
  edgeFile = paste("CytoscapeInput-edges-", paste(module, collapse="-"), ".txt", sep=""),
  nodeFile = paste("CytoscapeInput-nodes-", paste(module, collapse="-"), ".txt", sep=""),
  weighted = TRUE,
  threshold = 0.005,
  nodeNames = modProbes, 
  nodeAttr = moduleColors[inModule]
);