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

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

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

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

什么是WGCNA

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

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

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

下面我們來(lái)看一下WGCNA的官方流程：

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

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

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

sc-WGCNA

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

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

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

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

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

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

這一點(diǎn)高通量的單細(xì)胞轉(zhuǎn)錄組一般是不怕的，細(xì)胞數(shù)一般在幾千呀忧。但是這也是一個(gè)挑戰(zhàn)：A）數(shù)據(jù)緯度高了計(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分析之前也是每個(gè)cluster隨機(jī)選一部分細(xì)胞構(gòu)成Pseudocell（局部bulk的方法）泞辐。怪不得我用原始的count矩陣做WGCNA的結(jié)果這么差呢笔横。

傳統(tǒng)地，探針集或基因可以通過(guò)均值咐吼、絕對(duì)中位差(MAD)或方差進(jìn)行過(guò)濾吹缔，因?yàn)榈捅磉_(dá)或不變的基因通常代表噪聲。用均值表達(dá)還是方差過(guò)濾是否更好尚有爭(zhēng)議锯茄，兩者都有優(yōu)缺點(diǎn)厢塘。不建議通過(guò)差異分析過(guò)濾掉基因。我們知道單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)的一個(gè)特點(diǎn)就是緯度高撇吞，數(shù)據(jù)稀疏俗冻，除了要考慮細(xì)胞的特殊處理之外，我們還可以過(guò)濾基因牍颈，如只用高変基因（FindVariableFeatures()）。這里請(qǐng)注意琅关，也是不推薦只選用某一群的差異基因做的煮岁，因?yàn)槟骋蝗旱牟町惢颍呀?jīng)是一個(gè)明顯的模塊了涣易，這樣做很可能只得到很少的模塊画机。

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

WGCNA in action

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

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

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

載入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ù)很多泡躯，遠(yuǎn)大于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個(gè)細(xì)胞咕别，我們?cè)侔鸦蜻^(guò)濾一下。

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

最終我們的表達(dá)譜是這樣的：

 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

處理好表達(dá)譜写穴，下面的流程就常見了惰拱。

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

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

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

## 篩選中位絕對(duì)偏差前75%的基因翔冀，至少M(fèi)AD大于0.01
## 篩選后會(huì)降低運(yùn)算量，也會(huì)失去部分信息
## 也可不做篩選披泪，使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

## 檢測(cè)缺失值
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

注意一下磷箕，在篩選的時(shí)候光戈，有些基因是要保留的，哪些基因呢戳杀？就是那些你關(guān)注的基因艾少。

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

這要是幾千個(gè)細(xì)胞的話卡乾，簡(jiǎn)直一團(tuán)黑。

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ù)及其返回的對(duì)象缚够，最佳的beta值就是sft$powerEstimate

par(mfrow = c(1,2))
cex1 = 0.9
# 橫軸是Soft threshold (power)幔妨，縱軸是無(wú)標(biāo)度網(wǎng)絡(luò)的評(píng)估參數(shù)，數(shù)值越高谍椅，
# 網(wǎng)絡(luò)越符合無(wú)標(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ù)β雏吭，左圖縱軸代表對(duì)應(yīng)的網(wǎng)絡(luò)中l(wèi)og(k)與log(p(k))相關(guān)系數(shù)的平方锁施。相關(guān)系數(shù)的平方越高，說(shuō)明該網(wǎng)絡(luò)越逼近無(wú)網(wǎng)路尺度的分布杖们。右圖的縱軸代表對(duì)應(yīng)的基因模塊中所有基因鄰接函數(shù)的均值悉抵。最佳的beta值就是sft$powerEstimate，已經(jīng)被保存到變量了胀莹，不需要知道具體是什么基跑，后面的代碼都用這個(gè)即可，在本例子里面是5描焰。

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

# 無(wú)向網(wǎng)絡(luò)在power小于15或有向網(wǎng)絡(luò)power小于30內(nèi)，沒(méi)有一個(gè)power值可以使
# 無(wú)標(biāo)度網(wǎng)絡(luò)圖譜結(jié)構(gòu)R^2達(dá)到0.8篱竭，平均連接度較高如在100以上力图，可能是由于
# 部分樣品與其他樣品差別太大。這可能由批次效應(yīng)掺逼、樣品異質(zhì)性或?qū)嶒?yàn)條件對(duì)
# 表達(dá)影響太大等造成吃媒。可以通過(guò)繪制樣品聚類查看分組信息和有無(wú)異常樣品吕喘。
# 如果這確實(shí)是由有意義的生物變化引起的赘那，也可以使用下面的經(jīng)驗(yàn)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: 上一步計(jì)算的軟閾值
# maxBlockSize: 計(jì)算機(jī)能處理的最大模塊的基因數(shù)量 (默認(rèn)5000)氯质；
#  4G內(nèi)存電腦可處理8000-10000個(gè)募舟，16G內(nèi)存電腦可以處理2萬(wàn)個(gè)，32G內(nèi)存電腦可
#  以處理3萬(wàn)個(gè)
#  計(jì)算資源允許的情況下最好放在一個(gè)block里面闻察。
# corType: pearson or bicor
# numericLabels: 返回?cái)?shù)字而不是顏色作為模塊的名字拱礁，后面可以再轉(zhuǎn)換為顏色
# saveTOMs：最耗費(fèi)時(shí)間的計(jì)算，存儲(chǔ)起來(lái)辕漂，供后續(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)

計(jì)算過(guò)程：

 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...

我們肯定對(duì)這個(gè)返回的對(duì)象比較感興趣啦,，挨個(gè)查看一下钉嘹。

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 

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

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

我們這個(gè)至少不是全灰的啊褪那。

# module eigengene, 可以繪制線圖，作為每個(gè)模塊的基因表達(dá)趨勢(shì)的展示
MEs = net$MEs

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

# 根據(jù)基因間表達(dá)量進(jìn)行聚類所得到的各模塊間的相關(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)

計(jì)算每個(gè)模塊的特征向量基因，為某一特定模塊第一主成分基因E武学。

 計(jì)算每個(gè)模塊的特征向量基因祭往，為某一特定模塊第一主成分基因E。代表了該模塊內(nèi)基因表達(dá)的整體水平
MEList = moduleEigengenes(dataExpr, colors = dynamicColors)
MEs = MEList$eigengenes
# 計(jì)算根據(jù)模塊特征向量基因計(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) 

畫出指定模塊表達(dá)量的熱圖：

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

由于細(xì)胞名稱被我們改變了硼补，原來(lái)的對(duì)象我們要構(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ò)模塊與各種表型的相關(guān)系數(shù)，可以很清楚的挑選自己感興趣的模塊進(jìn)行下游分析了沛鸵。

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來(lái)挖掘與性狀相關(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ò)顯著性括勺，也就是性狀與每個(gè)基因表達(dá)量相關(guān)性在各個(gè)模塊的均值作為該性狀在該模塊的顯著性，顯著性最大的那個(gè)模塊與該性狀最相關(guān)：

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

# 獲得該性狀在每個(gè)模塊中的顯著性：
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）,首先計(jì)算基因的內(nèi)部連接度和模塊身份曲掰，內(nèi)部連接度衡量的是基因在模塊內(nèi)部的地位疾捍，而模塊身份表明基因?qū)儆谀膫€(gè)模塊。

# 計(jì)算每個(gè)基因模塊內(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)基于特征向量基因連接度評(píng)估模塊內(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)基于特征向量基因連接度評(píng)估模塊內(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]
);