作者：oct
審稿：童蒙
編輯：amethyst

一倒得、基本概念

WGCNA（ Weighted correlation network analysis）譯為加權(quán)基因共表達網(wǎng)絡分析第岖，是一種經(jīng)常用于基因組和系統(tǒng)生物學研究中的網(wǎng)絡分析方法馆蠕。在對一般的數(shù)據(jù)挖掘分享中沿彭，用來發(fā)現(xiàn)大量變量的成對關系瞻鹏，例如利用基因表達數(shù)據(jù)，來構(gòu)建基因共表達網(wǎng)絡捺萌，得到基因互作的關系档冬。

在許多WGCNA中，需要對網(wǎng)絡中各個模塊的特性桃纯，以及隨著不同情況下的變化酷誓，模塊的變化情況進行研究。比如态坦，可以研究不同組織或不同物種之間的模塊保存程度（module preservation）盐数。模塊保存程度統(tǒng)計信息（Preservation Statistic）旨在量化參考數(shù)據(jù)集和測試數(shù)據(jù)集之間模塊保存情況。WGCNA R包中的函數(shù)modulePreservation用于計算不同數(shù)據(jù)集之間的模塊保存統(tǒng)計信息伞梯。modulePreservation函數(shù)用法可以參考網(wǎng)頁：https://rdrr.io/cran/WGCNA/man/modulePreservation.html玫氢。

以往的模塊保存程度統(tǒng)計信息可以分為四大類：

• 交叉列表統(tǒng)計信息（Cross-tabulation）：會比較參考數(shù)據(jù)集和測試數(shù)據(jù)集中的聚類情況（即模塊類別分配），因此它們需要將聚類過程也應用于測試數(shù)據(jù)谜诫。
• 密度統(tǒng)計信息（Density）：不需要在測試數(shù)據(jù)集中分配模塊標簽漾峡。
• 可分離性統(tǒng)計信息（separability）：不需要在測試數(shù)據(jù)集中分配模塊標簽。
• 穩(wěn)定性統(tǒng)計信息（ stability）：當將一定數(shù)量的人工噪聲添加到數(shù)據(jù)中時喻旷，穩(wěn)定性統(tǒng)計通常研究數(shù)據(jù)集的穩(wěn)定性生逸，很少用在量化模塊保存程度。

在以往分類的基礎上掰邢，基于量化模塊保存度牺陶，目前又細分為以下五類：

• 交叉列表統(tǒng)計信息（Cross-tabulation）；
• 密度統(tǒng)計信息（Density）辣之；
• 可分離性統(tǒng)計信息（separability）掰伸；
• 連通性統(tǒng)計信息（Connectivity）：用于確定參考網(wǎng)絡基因之間的連接模式是否類似于測試網(wǎng)絡中的連接模式；
• 綜合統(tǒng)計信息（Composite）：整合密度統(tǒng)計信息與連通性統(tǒng)計信息怀估，綜合評估模塊保存度狮鸭。

這些統(tǒng)計數(shù)據(jù)通常通過每個數(shù)據(jù)集中相似對象（通常是基因）的數(shù)量或每個數(shù)據(jù)集之間不相似對象的數(shù)量來評估數(shù)據(jù)集。

modulePreservation函數(shù)多搀，一般會統(tǒng)計以下列表的統(tǒng)計信息歧蕉，其中Type列為每類統(tǒng)計量的類別，Network列為每類統(tǒng)計量使用的網(wǎng)絡類型康铭，*input列表示需要輸入的數(shù)據(jù)類型惯退，輸入數(shù)據(jù)需要以下類別：

• Lbl：module label；
• Adj：general network adjacency从藤；
• datX：numeric data from which a correlation network is constructed催跪。

二锁蠕、適用的場景

轉(zhuǎn)錄組或芯片數(shù)據(jù)在數(shù)據(jù)庫中存在大量類似的研究結(jié)果，那么這些結(jié)果之間的相容性如何懊蒸？WGCNA中modulePreservation函數(shù)用于闡述不同實驗結(jié)果的基因表達模式的關聯(lián)荣倾；

不同物種之間基因表達模式的分析也是多數(shù)文章的研究方向，WGCNA中modulePreservation函數(shù)同樣也適用于闡述不同物種的基因表達模式的關聯(lián)骑丸，比如人與黑猩猩基因表達模塊之間的基因保存性舌仍，類似的分析有很多，比如人與牛通危，人與豬等铸豁。

三、modulePreservation 函數(shù)實踐

以下命令用于演示在黑猩猩數(shù)據(jù)集中人腦共表達模塊的保存情況菊碟。

1. 配置R運行環(huán)境

啟動R之后推姻，設置工作目錄并加載所需的包。

# Display the current working directory
getwd()
# If necessary, change the path below to the directory where the data files are stored.
# "." means current directory. On Windows use a forward slash / instead of the usual 
workingDir = "."
setwd(workingDir)
# Load the package
library(WGCNA);
# The following setting is important, do not omit.
options(stringsAsFactors = FALSE)

2. 導入輸入數(shù)據(jù)

2.1 輸入表達矩陣

首先下載數(shù)據(jù)框沟，存放在本地目錄中藏古，表達矩陣數(shù)據(jù)存放在在Dataset 1 (network construction).csv.bz2中。

# data input
file = bzfile("Dataset 1 (network construction).csv.bz2")
dat1 = read.csv(file, header=T)

Dataset 1 (network construction).csv數(shù)據(jù)展示：

注：為展示方便忍燥，該數(shù)據(jù)共39列拧晕，中間省略32列。

表達矩陣包含36個樣本梅垄，每一行對應一個基因/探針厂捞，每一列對應一個樣本或輔助信息，我們需要刪除輔助信息队丝，并轉(zhuǎn)置表達矩陣用于后續(xù)分析靡馁。我們只保留BrainvariantH列>0的基因/探針。

dim(dat1)
names(dat1)
datExpr=data.frame(t(dat1[dat1$Brain_variant_H>0,2:39]))
indexHuman=c(19:36)
indexChimp=c(1:18)

現(xiàn)在我們設置多組表達數(shù)據(jù)和相應的模塊顏色机久。

# Number of data sets that we work with
nSets = 2
# Object that will contain the expression data
multiExpr = list()
multiExpr[[1]] = list(data = datExpr[indexHuman, ])
multiExpr[[2]] = list(data = datExpr[indexChimp, ])
# Names for the two sets
setLabels = c("Human", "Chimp")
# Important: components of multiExpr must carry identificating names
names(multiExpr) = setLabels
# Display the dimensions of the expression data (if you are confused by this construct, ignore it):
lapply(multiExpr, lapply, dim)

最后一條命令的輸出確認有兩個表達數(shù)據(jù)集臭墨，每個數(shù)據(jù)集包含18個樣本中4000個基因的表達測量值。

注意multiExpr的結(jié)構(gòu)：它是一個列表膘盖，每個輸入數(shù)據(jù)集必須包含一個稱為data的組件并且包含表達數(shù)據(jù)胧弛。外部列表必須具有適當?shù)拿Q，以便表達式數(shù)據(jù)可以與我們接下來創(chuàng)建的模塊標簽匹配侠畔。

2.2 加載模塊標簽

WGCNA分析中结缚，通過構(gòu)建共表達矩陣，計算基因間的相似性软棺，確定基因模塊红竭，通常用不同的顏色來代表模塊，即每個模塊中的基因?qū)环Ncolor，方便后續(xù)的模塊可視化茵宪。后續(xù)會分享該部分的相關內(nèi)容斥黑。

HumanChimp-OldhamAnalysis-colorHuman-colorChimp-inNetwork.RData中包含黑猩猩表達矩陣中每個基因的模塊顏色信息（colorChimp）以及人表達矩陣中每個基因的模塊顏色信息（colorHuman）。

# 加載從網(wǎng)絡分析獲得的模塊標簽眉厨，并創(chuàng)建一個列表。
x = load("HumanChimp-OldhamAnalysis-colorHuman-colorChimp-inNetwork.RData") # Create an object (list) holding the module labels for each set:
colorList = list(colorHuman, colorChimp)
# Components of the list must be named so that the names can be matched to the names of multiExpr
names(colorList) = setLabel

注意兽狭，multiExpr的名稱和multiColor（即這里的colorList）的名稱必須對應憾股。

3. 對模塊保存程度進行統(tǒng)計分析

這里使用modulePreservation進行模塊保存性計算，運行后將結(jié)果保存為Rdata箕慧，下次可直接加載服球，節(jié)省時間。

system.time( {
 mp = modulePreservation(multiExpr, colorList,
             referenceNetworks = c(1:2),
             loadPermutedStatistics = FALSE,
             nPermutations = 200,
             verbose = 3)
} )

save(mp, file = "HumanChimp-HumanSpecific-modulePreservation.RData")

該函數(shù)輸出的是一個嵌套列表颠焦，包括quality斩熊，preservation，accracy伐庭，permutationDetails粉渠，referenceSeparability和testSeparability，每個列表分別包含4或5個組成部分圾另，包含：

• observed（觀察值）霸株；
• Z（Z得分）；
• log.p（以10為底的p值的對數(shù)）集乔；
• log.pBonf（Bonferoni校正后的p值的以10為底的對數(shù)）去件；
• q（可選，q值）扰路。

每個列表包含Z尤溜，log.p，log.pBonf（可選的q值）汗唱，accuracy中包含的observedOverlapCounts和observedFisherPvalues被構(gòu)造為2級列表宫莱，其中外部組件對應于參考集，內(nèi)部組件對應于測試集哩罪。例如梢睛，preservation$observed[[1]] [[2]]包含用于保存測試集中參考集的密度和保存性等保存統(tǒng)計信息，1是參考集识椰，2是測試集绝葡。

4. 對模塊保存程度進行結(jié)果展示

ref = 1 # Select the human data as reference
test = 2 # Select the chimp data as test
statsObs = cbind(mp$quality$observed[[ref]][[test]][, -1], mp$preservation$observed[[ref]][[test]][, -1])
statsZ = cbind(mp$quality$Z[[ref]][[test]][, -1], mp$preservation$Z[[ref]][[test]][, -1])

我們看主要輸出Zsummary評分。

print(signif(statsZ[, "Zsummary.pres", drop = FALSE],2))

一般文獻中認為：

• Zsummary>10：高度保存腹鹉；
• 2<Zsummary<10：弱保存藏畅；
• Zsummary<2：無保存證據(jù)。

查看不同物種的不同模塊之間的overlap基因數(shù)與overlap pvalue。

overlap = cbind(mp$accuracy$observedCounts[[1]][[2]], mp$accuracy$observedFisherPvalues[[1]][[2]])
print(overlap)

輸出結(jié)果的前8列為overlap基因數(shù)愉阎，后8列為 overlap pvalue绞蹦。

參考文獻

• Jiang Z, Sun J, Dong H, Luo O, Zheng X, Obergfell C, Tang Y, Bi J, O'Neill R, Ruan Y, Chen J, Tian XC. Transcriptional profiles of bovine in vivo pre-implantation development. BMC Genomics. 2014 Sep 4;15(1):756. doi: 10.1186/1471-2164-15-756. PMID: 25185836; PMCID: PMC4162962.
• Langfelder P, Luo R, Oldham MC, Horvath S. Is my network module preserved and reproducible? PLoS Comput Biol. 2011 Jan 20;7(1):e1001057. doi: 10.1371/journal.pcbi.1001057. PMID: 21283776; PMCID: PMC3024255.
• Khan FA, Liu H, Zhou H, Wang K, Qamar MTU, Pandupuspitasari NS, Shujun Z. Analysis of Bos taurus and Sus scrofa X and Y chromosome transcriptome highlights reproductive driver genes. Oncotarget. 2017 Apr 13;8(33):54416-54433. doi: 10.18632/oncotarget.17081. PMID: 28903352; PMCID: PMC5589591.
• Li Y, Sun J, Ling Y, Ming H, Chen Z, Fang F, Liu Y, Cao H, Ding J, Cao Z, Zhang X, Bondioli K, Jiang Z, Zhang Y. Transcription profiles of oocytes during maturation and embryos during preimplantation development in vivo in the goat. Reprod Fertil Dev. 2020 Apr;32(7):714-725. doi: 10.1071/RD19391. PMID: 32317096.
• Xue Z, Huang K, Cai C, Cai L, Jiang CY, Feng Y, Liu Z, Zeng Q, Cheng L, Sun YE, Liu JY, Horvath S, Fan G. Genetic programs in human and mouse early embryos revealed by single-cell RNA sequencing. Nature. 2013 Aug 29;500(7464):593-7. doi: 10.1038/nature12364. Epub 2013 Jul 28. PMID: 23892778; PMCID: PMC4950944.
• Miller JA, Horvath S, Geschwind DH. Divergence of human and mouse brain transcriptome highlights Alzheimer disease pathways. Proc Natl Acad Sci U S A. 2010 Jul * 13;107(28):12698-703. doi: 10.1073/pnas.0914257107. Epub 2010 Jun 25. PMID: 20616000; PMCID: PMC2906579.