WGCNA怎么挖掘潛在的共表達基因

這篇文章主要介紹“WGCNA怎么挖掘潛在的共表達基因”，在日常操作中，相信很多人在WGCNA怎么挖掘潛在的共表達基因問題上存在疑惑，小編查閱了各式資料，整理出簡單好用的操作方法，希望對大家解答”WGCNA怎么挖掘潛在的共表達基因”的疑惑有所幫助！接下來，請跟著小編一起來學習吧！

共表達基因指的是表達量具有協同變化趨勢的基因集合，通常認為這些基因參與相同的生物學過程，比如參與同一個代謝通路，正是由于功能上的協同作用，導致表達量呈現出高度相關性。

在WGCNA中，對傳統的相關系數進行乘方運算，用最終得到的值來表征基因間的相關性。在計算出這樣的相關性統計量值之后，如何確定哪些基因是共表達的呢？

WGCNA的做法是聚類分析，聚類分析屬于一種非監督的機器學習算法，通過聚類樹，可以觀察到哪些基因在聚類樹中屬于同一分支，屬于同一分支的基因可以歸為一類。實際操作中，考慮到基因數目較多等情況，肯定需要算法來自動化的進行分類，WGCNA采用的是dynamicTreeCut這個R包。

對于聚類算法而言，需要輸入基因間的距離矩陣，首先就需要將基因間的鄰接矩陣轉換為距離矩陣，對相關系數進行乘方運算，可以計算出鄰接矩陣，但是這個值本質上反映的是基因間的相似度，并不是距離。在計算距離矩陣時，WGCNA采用了TOM這種統計量，該統計量可以表征網絡中節點的相似性，計算公式如下

對于兩個基因i和j而言，a表示兩個基因鄰接矩陣中對應的值，就是相關系數的乘方，K代表的每個基因的連接度， 公式如下

對于加權網絡而言，就是該節點的邊對應數值的總和，比如在網絡中基因A與3個基因相連，基因A的連接度就是對應3條邊的數值之和。兩個基因間的l值代表的是兩個基因所有邊的權重乘積的總和，公式如下

公式只是幫助我們理解計算的過程，其實只需要理解TOM是表征節點的相似度就行，我們要的是距離，所以直接用1減去相似度即可，公式如下

借助TOM值，將基因間的相關系數轉換為了距離，然后就可以用該距離矩陣進行聚類。上述的計算方法在WGCNA中都有對應的公式，代碼如下

# 確定乘方運算中power的最佳取值
powers <- c(c(1:10), seq(from = 12, to=20, by=2))

sft <- pickSoftThreshold(
datExpr,
powerVector = powers,
verbose = 5)

softPower <- sft$powerEstimate

#  計算鄰接矩陣
adjacency <- adjacency(datExpr, power = softPower)

# 計算TOM相似度矩陣
TOM <- TOMsimilarity(adjacency)

# 計算距離矩陣
dissTOM <- 1-TOM

# 聚類
geneTree <- hclust(as.dist(dissTOM), method = "average")

根據聚類結果和距離矩陣，就可以調用dynamicTreeCut的算法來識別modules, 代碼如下

# 指定每個module中基因數目的最小值
minModuleSize <- 30

# 識別modules
dynamicMods <- cutreeDynamic(
dendro = geneTree,
distM = dissTOM,
deepSplit = 2,
pamRespectsDendro = FALSE,
minClusterSize = minModuleSize)

通過table函數可以查看modules的結果，用法如下

> table(dynamicMods)
dynamicMods
  0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19
88 614 316 311 257 235 225 212 158 153 121 106 102 100  94  91  78  76  65  58
20  21  22
58  48  34

可以看到，識別出22個modules, 0代表那些沒有歸入任何modules的基因。通過plotDendroAndColors函數可視化聚類樹對應的modules, 代碼如下

dynamicColors = labels2colors(dynamicMods)
plotDendroAndColors(
geneTree,
dynamicColors,
"Dynamic Tree Cut",
dendroLabels = FALSE, hang = 0.03,
addGuide = TRUE,
guideHang = 0.05,
main = "Gene dendrogram and module colors")

生成的圖片如下

整個圖片分為兩個部分，上方為基因的聚類樹，下方為識別到的modules， 不同的modules對應不同的顏色，其中灰色對應那些沒有歸入任何modules的基因。

通過dynamicTreeCut識別到modules之后，還會結合每個modules的基因表達量數據，來識別相關性很高的modules, 從而進行合并，其原理是對modules進行聚類，每個module下的基因表達量是一個二維矩陣，做相關性分析我們只需要一個一維向量就可以了，可以利用PCA分析提取第一主成分來表征原始的矩陣，在WGCNA中，把每個module的表達譜數據對應的一維向量稱之為Module eigengene E。獲取一維向量之后，就可以計算相關性，直接用1減去相關性作為距離，進行聚類，代碼如下

MEList <- moduleEigengenes(
  datExpr, 
   colors = dynamicColors)
MEs <- MEList$eigengenes
MEDiss <- 1-cor(MEs)
METree <- hclust(as.dist(MEDiss), method = "average")
plot(METree,
main = "Clustering of module eigengenes",
xlab = "", sub = "")

modules的聚類樹示意如下

每個modules的名字用對應的顏色表示，在該聚類數中，分支長度為1減去兩個module間的相關系數，在合并modules時，將高相關性的合并為一類，可以指定一個閾值，比如將相關系數大于0.8的合并為一類，在該聚類樹中，對應的就是height小于0.2的modules,  對應下圖紅色的線

可以看到有8個modules都滿足條件，在合并時，會將原本屬于同一分支的modules直接合并為一個，從圖上可以看出，合并后會減少4個modules。合并的代碼如下

MEDissThres <- 0.2
merge <- mergeCloseModules(
  datExpr,
  dynamicColors,
  cutHeight = MEDissThres,
  verbose = 3)
mergedColors <- merge$colors
mergedMEs <- merge$newMEs
plotDendroAndColors(
geneTree,
cbind(dynamicColors, mergedColors),
c("Dynamic Tree Cut", "Merged dynamic"),
dendroLabels = FALSE,
hang = 0.03,
addGuide = TRUE,
guideHang = 0.05)

合并之后的modules 對應的圖片如下

到此，關于“WGCNA怎么挖掘潛在的共表達基因”的學習就結束了，希望能夠解決大家的疑惑。理論與實踐的搭配能更好的幫助大家學習，快去試試吧！若想繼續學習更多相關知識，請繼續關注億速云網站，小編會繼續努力為大家帶來更多實用的文章！