WGCNA 分析
基本概念
WGCNA其譯為加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析。該分析方法旨在尋找協(xié)同表達(dá)的基因模塊(module)，并探索基因網(wǎng)絡(luò)與關(guān)注的表型之間的關(guān)聯(lián)關(guān)系，以及網(wǎng)絡(luò)中的核心基因。
適用于復(fù)雜的數(shù)據(jù)模式，推薦5組(或者15個(gè)樣品)以上的數(shù)據(jù)。一般可應(yīng)用的研究方向有：不同器官或組織類型發(fā)育調(diào)控、同一組織不同發(fā)育調(diào)控、非生物脅迫不同時(shí)間點(diǎn)應(yīng)答、病原菌侵染后不同時(shí)間點(diǎn)應(yīng)答。
基本原理
從方法上來講，WGCNA分為表達(dá)量聚類分析和表型關(guān)聯(lián)兩部分，主要包括基因之間相關(guān)系數(shù)計(jì)算、基因模塊的確定、共表達(dá)網(wǎng)絡(luò)、模塊與性狀關(guān)聯(lián)四個(gè)步驟。
第一步計(jì)算任意兩個(gè)基因之間的相關(guān)系數(shù)（Person Coefficient）。為了衡量兩個(gè)基因是否具有相似表達(dá)模式，一般需要設(shè)置閾值來篩選，高于閾值的則認(rèn)為是相似的。但是這樣如果將閾值設(shè)為0.8，那么很難說明0.8和0.79兩個(gè)是有顯著差別的。因此，WGCNA分析時(shí)采用相關(guān)系數(shù)加權(quán)值，即對(duì)基因相關(guān)系數(shù)取N次冪，使得網(wǎng)絡(luò)中的基因之間的連接服從無尺度網(wǎng)絡(luò)分布(scale-freenetworks)，這種算法更具生物學(xué)意義。
第二步通過基因之間的相關(guān)系數(shù)構(gòu)建分層聚類樹，聚類樹的不同分支代表不同的基因模塊，不同顏色代表不同的模塊。基于基因的加權(quán)相關(guān)系數(shù)，將基因按照表達(dá)模式進(jìn)行分類，將模式相似的基因歸為一個(gè)模塊。這樣就可以將幾萬個(gè)基因通過基因表達(dá)模式被分成了幾十個(gè)模塊，是一個(gè)提取歸納信息的過程。
WGCNA術(shù)語
權(quán)重(weghted)
基因之間不僅僅是相關(guān)與否，還記錄著它們的相關(guān)性數(shù)值，數(shù)值就是基因之間的聯(lián)系的權(quán)重(相關(guān)性)。
[圖片上傳中...(image-ae7c41-1616404287726-8)]
Module
模塊(module)：表達(dá)模式相似的基因分為一類，這樣的一類基因成為模塊；
Eigengene
Eigengene（eigen + gene）：基因和樣本構(gòu)成的矩陣，https://en.wiktionary.org/wiki/eigengene
Adjacency Matrix
鄰近矩陣：是圖的一種存儲(chǔ)形式，用一個(gè)一維數(shù)組存放圖中所有頂點(diǎn)數(shù)據(jù)；用一個(gè)二維數(shù)組存放頂點(diǎn)間關(guān)系（邊或?。┑臄?shù)據(jù)，這個(gè)二維數(shù)組稱為鄰接矩陣；在WGCNA分析里面指的是基因與基因之間的相關(guān)性系數(shù)矩陣。 如果用了閾值來判斷基因相關(guān)與否，那么這個(gè)鄰近矩陣就是0/1矩陣，只記錄基因相關(guān)與否。但是WGCNA沒有用閾值來卡基因的相關(guān)性，而是記錄了所有基因之間的相關(guān)性。
Topological Overlap Matrix (TOM)
WGNA認(rèn)為基因之間的簡單的相關(guān)性不足以計(jì)算共表達(dá)，所以它利用上面的鄰近矩陣，又計(jì)算了一個(gè)新的鄰近矩陣。一般來說，TOM就是WGCNA分析的最終結(jié)果，后續(xù)的只是對(duì)TOM的下游注釋。
下游分析
得到模塊之后的分析有：
1.模塊的功能富集
2.模塊與性狀之間的相關(guān)性
3.模塊與樣本間的相關(guān)系數(shù)
挖掘模塊的關(guān)鍵信息：
1.找到模塊的核心基因
2.利用關(guān)系預(yù)測(cè)基因功能
代碼示例
其中第一步數(shù)據(jù)準(zhǔn)備反而是最復(fù)雜的，取決于大家的R語言水平，這個(gè)數(shù)據(jù)GSE48213-wgcna-input.RData我已經(jīng)保存下來咯，如果大家不會(huì)做，又想體驗(yàn)一下這個(gè)WGCNA流程，就可以直接load我保存好的數(shù)據(jù)文件即可。
step1: 輸入數(shù)據(jù)的準(zhǔn)備
這里主要是表達(dá)矩陣， 如果是芯片數(shù)據(jù)，那么常規(guī)的歸一化矩陣即可，如果是轉(zhuǎn)錄組數(shù)據(jù)，最好是RPKM/TPM值或者其它歸一化好的表達(dá)量。然后就是臨床信息或者其它表型，總之就是樣本的屬性。
為了保證后續(xù)腳本的統(tǒng)一性，表達(dá)矩陣統(tǒng)一用datExpr標(biāo)識(shí)，臨床 信息統(tǒng)一用datTraits標(biāo)識(shí)。(PS: 如果你R語言很差，變量名不要輕易修改)
library(WGCNA)
RNAseq_voom <- fpkm
## 因?yàn)閃GCNA針對(duì)的是基因進(jìn)行聚類，而一般我們的聚類是針對(duì)樣本用hclust即可，所以這個(gè)時(shí)候需要轉(zhuǎn)置。
WGCNA_matrix = t(RNAseq_voom[order(apply(RNAseq_voom,1,mad), decreasing = T)[1:5000],])
datExpr0 <- WGCNA_matrix ## top 5000 mad genes
datExpr <- datExpr0
## 下面主要是為了防止臨床表型與樣本名字對(duì)不上
sampleNames = rownames(datExpr);
traitRows = match(sampleNames, datTraits$gsm)
rownames(datTraits) = datTraits[traitRows, 1]
上面代碼里面的rpkm就是我們的轉(zhuǎn)錄組數(shù)據(jù)的表達(dá)矩陣，以rpkm為單位。而datTraits就是所有樣本對(duì)應(yīng)的表型信息。需要自己制作，這個(gè)是學(xué)習(xí)WGCNA的基礎(chǔ)，本次實(shí)例代碼都是基于這兩個(gè)數(shù)據(jù)。至于如何做出上面代碼的兩個(gè)例子，取決于大家自己的項(xiàng)目，我這里給出自己的代碼，僅供參考哈！
setwd('WGCNA/')
# 56 breast cancer cell lines were profiled to identify patterns of gene expression associated with subtype and response to therapeutic compounds.
if(F){
## https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE48213
#wget -c ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE48nnn/GSE48213/suppl/GSE48213_RAW.tar
#tar -xf GSE48213_RAW.tar
#gzip -d *.gz
## 首先在GSE48213_RAW目錄里面生成tmp.txt文件，使用shell腳本
# awk '{print FILENAME"\t"$0}' GSM*.txt |grep -v EnsEMBL_Gene_ID >tmp.txt
# 其實(shí)也可以直接使用R來讀取GSE48213_RAW.tar里面的gz文件，這里就不演示了
# 可以參考：https://mp.weixin.qq.com/s/OLc9QmfN0YcT548VAYgOPA 里面的教程
## 然后把tmp.txt導(dǎo)入R語言里面用reshape2處理即可
# 這個(gè) tmp.txt 文件應(yīng)該是100M左右大小哦。
a=read.table('GSE48213_RAW/tmp.txt',sep = '\t',stringsAsFactors = F)
library(reshape2)
fpkm <- dcast(a,formula = V2~V1)
rownames(fpkm)=fpkm[,1]
fpkm=fpkm[,-1]
colnames(fpkm)=sapply(colnames(fpkm),function(x) strsplit(x,"_")[[1]][1])
library(GEOquery)
a=getGEO('GSE48213')
metadata=pData(a[[1]])[,c(2,10,12)]
datTraits = data.frame(gsm=metadata[,1],
cellline=trimws(sapply(as.character(metadata$characteristics_ch1),function(x) strsplit(x,":")[[1]][2])),
subtype=trimws(sapply(as.character(metadata$characteristics_ch1.2),function(x) strsplit(x,":")[[1]][2]))
)
save(fpkm,datTraits,file = 'GSE48213-wgcna-input.RData')
}
很明顯，這個(gè)數(shù)據(jù)GSE48213-wgcna-input.RData我已經(jīng)保存下來咯，如果大家不會(huì)做，又想體驗(yàn)一下這個(gè)WGCNA流程，那么可以找我email求取這個(gè)數(shù)據(jù)哦。我的郵箱是jmzeng1314@163.com
我給大家演示的示例數(shù)據(jù)大概是下面這個(gè)樣子：
> head(datTraits) ## 56 個(gè)細(xì)胞系的分類信息，表型
gsm cellline subtype
GSM1172844 GSM1172844 184A1 Non-malignant
GSM1172845 GSM1172845 184B5 Non-malignant
GSM1172846 GSM1172846 21MT1 Basal
GSM1172847 GSM1172847 21MT2 Basal
GSM1172848 GSM1172848 21NT Basal
GSM1172849 GSM1172849 21PT Basal
> fpkm[1:4,1:4] ## 56個(gè)細(xì)胞系的36953個(gè)基因的表達(dá)矩陣
GSM1172844 GSM1172845 GSM1172846 GSM1172847
ENSG00000000003 95.21255 95.69868 19.99467 65.6863763
ENSG00000000005 0.00000 0.00000 0.00000 0.1492021
ENSG00000000419 453.20831 243.64804 142.05818 200.4131493
ENSG00000000457 18.10439 26.56661 16.12776 12.0873135
>
這個(gè)數(shù)據(jù)集里面的56種細(xì)胞系被分成了5組，如果要分開兩兩做差異分析，有10種組合，也就是說需要做10次差異分析，每個(gè)差異分析結(jié)果都需要去注釋，會(huì)比較麻煩，這個(gè)時(shí)候WGCNA就派上用場啦。當(dāng)然，如果你一定要去做差異分析，我也給你代碼：https://github.com/jmzeng1314/my-R/blob/master/10-RNA-seq-3-groups/hisat2_mm10_htseq.R
實(shí)際上多個(gè)分組，差異分析策略是非常個(gè)性化的， 比如：https://mp.weixin.qq.com/s/hc6JkKxyelc7b1M1MRiHRQ
step2:確定最佳beta值
選擇合適“軟閥值（soft thresholding power）”beta，同樣的，也是使用教程標(biāo)準(zhǔn)代碼即可：
powers = c(c(1:10), seq(from = 12, to=20, by=2))
# Call the network topology analysis function
sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)
#設(shè)置網(wǎng)絡(luò)構(gòu)建參數(shù)選擇范圍，計(jì)算無尺度分布拓?fù)渚仃?
# Plot the results:
##sizeGrWindow(9, 5)
par(mfrow = c(1,2));
cex1 = 0.9;
# Scale-free topology fit index as a function of the soft-thresholding power
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")
# 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")
關(guān)鍵就是理解pickSoftThreshold函數(shù)及其返回的對(duì)象，最佳的beta值就是sft$powerEstimate
參數(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ù)的平方越高，說明該網(wǎng)絡(luò)越逼近無網(wǎng)路尺度的分布。右圖的縱軸代表對(duì)應(yīng)的基因模塊中所有基因鄰接函數(shù)的均值。最佳的beta值就是sft$powerEstimate，已經(jīng)被保存到變量了，不需要知道具體是什么，后面的代碼都用這個(gè)即可，在本例子里面是6。
即使你不理解它，也可以使用代碼拿到合適“軟閥值（soft thresholding power）”beta進(jìn)行后續(xù)分析。
step3：一步法構(gòu)建共表達(dá)矩陣
有了表達(dá)矩陣和估計(jì)好的最佳beta值，就可以直接構(gòu)建共表達(dá)矩陣了。
net = blockwiseModules(
datExpr,
power = sft$powerEstimate,
maxBlockSize = 6000,
TOMType = "unsigned", minModuleSize = 30,
reassignThreshold = 0, mergeCutHeight = 0.25,
numericLabels = TRUE, pamRespectsDendro = FALSE,
saveTOMs = F,
verbose = 3
)
table(net$colors)
所有的核心就在這一步，把輸入的表達(dá)矩陣的幾千個(gè)基因組歸類成了幾十個(gè)模塊。大體思路：計(jì)算基因間的鄰接性，根據(jù)鄰接性計(jì)算基因間的相似性，然后推出基因間的相異性系數(shù)，并據(jù)此得到基因間的系統(tǒng)聚類樹。然后按照混合動(dòng)態(tài)剪切樹的標(biāo)準(zhǔn)，設(shè)置每個(gè)基因模塊最少的基因數(shù)目為30。
根據(jù)動(dòng)態(tài)剪切法確定基因模塊后，再次分析，依次計(jì)算每個(gè)模塊的特征向量值，然后對(duì)模塊進(jìn)行聚類分析，將距離較近的模塊合并為新的模塊。
step4: 模塊可視化
這里用不同的顏色來代表那些所有的模塊，其中灰色默認(rèn)是無法歸類于任何模塊的那些基因，如果灰色模塊里面的基因太多，那么前期對(duì)表達(dá)矩陣挑選基因的步驟可能就不太合適。
# Convert labels to colors for plotting
mergedColors = labels2colors(net$colors)
table(mergedColors)
# Plot the dendrogram and the module colors underneath
plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
"Module colors",
dendroLabels = FALSE, hang = 0.03,
addGuide = TRUE, guideHang = 0.05)
## assign all of the gene to their corresponding module
## hclust for the genes.
這里的重點(diǎn)就是plotDendroAndColors函數(shù)，它接受一個(gè)聚類的對(duì)象，以及該對(duì)象里面包含的所有個(gè)體所對(duì)應(yīng)的顏色。比如對(duì)表達(dá)矩陣進(jìn)行hclust之后，加上表達(dá)矩陣?yán)锩嫠袠颖镜姆纸M信息對(duì)應(yīng)的顏色，也是可以用plotDendroAndColors函數(shù)可視化的，比如下面的樣品圖：
#明確樣本數(shù)和基因數(shù)
nGenes = ncol(datExpr)
nSamples = nrow(datExpr)
#首先針對(duì)樣本做個(gè)系統(tǒng)聚類樹
datExpr_tree<-hclust(dist(datExpr), method = "average")
par(mar = c(0,5,2,0))
plot(datExpr_tree, main = "Sample clustering", sub="", xlab="", cex.lab = 2,
cex.axis = 1, cex.main = 1,cex.lab=1)
## 如果這個(gè)時(shí)候樣本是有性狀，或者臨床表型的，可以加進(jìn)去看看是否聚類合理
#針對(duì)前面構(gòu)造的樣品矩陣添加對(duì)應(yīng)顏色
sample_colors <- numbers2colors(as.numeric(factor(datTraits$Tumor.Type)),
colors = c("white","blue","red","green"),signed = FALSE)
## 這個(gè)給樣品添加對(duì)應(yīng)顏色的代碼需要自行修改以適應(yīng)自己的數(shù)據(jù)分析項(xiàng)目。
# sample_colors <- numbers2colors( datTraits ,signed = FALSE)
## 如果樣品有多種分類情況，而且 datTraits 里面都是分類信息，那么可以直接用上面代碼，當(dāng)然，這樣給的顏色不明顯，意義不大。
#構(gòu)造10個(gè)樣品的系統(tǒng)聚類樹及性狀熱圖
par(mar = c(1,4,3,1),cex=0.8)
plotDendroAndColors(datExpr_tree, sample_colors,
groupLabels = colnames(sample),
cex.dendroLabels = 0.8,
marAll = c(1, 4, 3, 1),
cex.rowText = 0.01,
main = "Sample dendrogram and trait heatmap")
上面給樣本進(jìn)行聚類的代碼可以不運(yùn)行，其實(shí)跟WGCNA本身關(guān)系不大。
可以看到這些乳腺癌的細(xì)胞系的表達(dá)譜聚類情況并不是完全與其分類匹配，所以僅僅是根據(jù)樣本的分組信息做差異分析并不完全準(zhǔn)確。
step5:模塊和性狀的關(guān)系
## step 5 (最重要的) 模塊和性狀的關(guān)系
## 這一步主要是針對(duì)于連續(xù)變量，如果是分類變量，需要轉(zhuǎn)換成連續(xù)變量方可使用
table(datTraits$subtype)
if(T){
nGenes = ncol(datExpr)
nSamples = nrow(datExpr)
design=model.matrix(~0+ datTraits$subtype)
colnames(design)=levels(datTraits$subtype)
moduleColors <- labels2colors(net$colors)
# Recalculate MEs with color labels
MEs0 = moduleEigengenes(datExpr, moduleColors)$eigengenes
MEs = orderMEs(MEs0); ##不同顏色的模塊的ME值矩 (樣本vs模塊)
moduleTraitCor = cor(MEs, design , use = "p");
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)
sizeGrWindow(10,6)
# Will display correlations and their p-values
textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
signif(moduleTraitPvalue, 1), ")", sep = "");
dim(textMatrix) = dim(moduleTraitCor)
png("step5-Module-trait-relationships.png",width = 800,height = 1200,res = 120)
par(mar = c(6, 8.5, 3, 3));
# Display the correlation values within a heatmap plot
labeledHeatmap(Matrix = moduleTraitCor,
xLabels = colnames(design),
yLabels = names(MEs),
ySymbols = names(MEs),
colorLabels = FALSE,
colors = greenWhiteRed(50),
textMatrix = textMatrix,
setStdMargins = FALSE,
cex.text = 0.5,
zlim = c(-1,1),
main = paste("Module-trait relationships"))
dev.off()
# 除了上面的熱圖展現(xiàn)形狀與基因模塊的相關(guān)性外
# 還可以是條形圖,但是只能是指定某個(gè)形狀
# 或者自己循環(huán)一下批量出圖。
Luminal = as.data.frame(design[,3]);
names(Luminal) = "Luminal"
y=Luminal
GS1=as.numeric(cor(y,datExpr, use="p"))
GeneSignificance=abs(GS1)
# Next module significance is defined as average gene significance.
ModuleSignificance=tapply(GeneSignificance,
moduleColors, mean, na.rm=T)
sizeGrWindow(8,7)
par(mfrow = c(1,1))
# 如果模塊太多，下面的展示就不友好
# 不過，我們可以自定義出圖。
plotModuleSignificance(GeneSignificance,moduleColors)
}
通過模塊與各種表型的相關(guān)系數(shù)，可以很清楚的挑選自己感興趣的模塊進(jìn)行下游分析了。這個(gè)圖就是把moduleTraitCor這個(gè)矩陣給用熱圖可視化一下。
因?yàn)橐恍v史遺留問題，這個(gè)圖片缺乏X軸的標(biāo)記。
從上圖已經(jīng)可以看到跟乳腺癌分類相關(guān)的基因模塊了，包括"Basal" "Claudin-low" "Luminal" "Non-malignant" "unknown"這5類所對(duì)應(yīng)的不同模塊的基因列表?？梢钥吹矫恳环N乳腺癌都有跟它強(qiáng)烈相關(guān)的模塊，可以作為它的表達(dá)signature，模塊里面的基因可以拿去做下游分析。我們看到Luminal表型跟棕色的模塊相關(guān)性高達(dá)0.86，而且極其顯著的相關(guān)，所以值得我們挖掘，這個(gè)模塊里面的基因是什么，為什么如此的相關(guān)呢？
step6:感興趣性狀的模塊的具體基因分析
性狀跟模塊雖然求出了相關(guān)性，可以挑選最相關(guān)的那些模塊來分析，但是模塊本身仍然包含非常多的基因，還需進(jìn)一步的尋找最重要的基因。所有的模塊都可以跟基因算出相關(guān)系數(shù)，所有的連續(xù)型性狀也可以跟基因的表達(dá)值算出相關(guān)系數(shù)。主要參考資料：PDF document, R script 如果跟性狀顯著相關(guān)基因也跟某個(gè)模塊顯著相關(guān)，那么這些基因可能就非常重要。
首先計(jì)算模塊與基因的相關(guān)性矩陣
# names (colors) of the modules
modNames = substring(names(MEs), 3)
geneModuleMembership = as.data.frame(cor(datExpr, MEs, use = "p"));
## 算出每個(gè)模塊跟基因的皮爾森相關(guān)系數(shù)矩陣
## MEs是每個(gè)模塊在每個(gè)樣本里面的值
## datExpr是每個(gè)基因在每個(gè)樣本的表達(dá)量
MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples));
names(geneModuleMembership) = paste("MM", modNames, sep="");
names(MMPvalue) = paste("p.MM", modNames, sep="");
再計(jì)算性狀與基因的相關(guān)性矩陣
## 只有連續(xù)型性狀才能只有計(jì)算
## 這里把是否屬于 Luminal 表型這個(gè)變量用0,1進(jìn)行數(shù)值化。
Luminal = as.data.frame(design[,3]);
names(Luminal) = "Luminal"
geneTraitSignificance = as.data.frame(cor(datExpr, Luminal, use = "p"));
GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples));
names(geneTraitSignificance) = paste("GS.", names(Luminal), sep="");
names(GSPvalue) = paste("p.GS.", names(Luminal), sep="");
最后把兩個(gè)相關(guān)性矩陣聯(lián)合起來,指定感興趣模塊進(jìn)行分析
module = "brown"
column = match(module, modNames);
moduleGenes = moduleColors==module;
sizeGrWindow(7, 7);
par(mfrow = c(1,1));
verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
abs(geneTraitSignificance[moduleGenes, 1]),
xlab = paste("Module Membership in", module, "module"),
ylab = "Gene significance for Luminal",
main = paste("Module membership vs. gene significance\n"),
cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
可以看到這些基因不僅僅是跟其對(duì)應(yīng)的模塊高度相關(guān)，而且是跟其對(duì)應(yīng)的性狀高度相關(guān)，進(jìn)一步說明了基因值得深度探究。
step7:網(wǎng)絡(luò)的可視化
主要參考資料：PDF document, R script
首先針對(duì)所有基因畫熱圖
# 主要是可視化 TOM矩陣，WGCNA的標(biāo)準(zhǔn)配圖
# 然后可視化不同 模塊 的相關(guān)性 熱圖
# 不同模塊的層次聚類圖
# 還有模塊診斷，主要是 intramodular connectivity
if(T){
nGenes = ncol(datExpr)
nSamples = nrow(datExpr)
geneTree = net$dendrograms[[1]];
dissTOM = 1-TOMsimilarityFromExpr(datExpr, power = 6);
plotTOM = dissTOM^7;
diag(plotTOM) = NA;
#TOMplot(plotTOM, geneTree, moduleColors, main = "Network heatmap plot, all genes")
nSelect = 400
# For reproducibility, we set the random seed
set.seed(10);
select = sample(nGenes, size = nSelect);
selectTOM = dissTOM[select, select];
# There’s no simple way of restricting a clustering tree to a subset of genes, so we must re-cluster.
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;
png("step7-Network-heatmap.png",width = 800,height = 600)
TOMplot(plotDiss, selectTree, selectColors, main = "Network heatmap plot, selected genes")
dev.off()
# Recalculate module eigengenes
MEs = moduleEigengenes(datExpr, moduleColors)$eigengenes
## 只有連續(xù)型性狀才能只有計(jì)算
## 這里把是否屬 Luminal 表型這個(gè)變量0,1進(jìn)行數(shù)值化
Luminal = as.data.frame(design[,3]);
names(Luminal) = "Luminal"
# Add the weight to existing module eigengenes
MET = orderMEs(cbind(MEs, Luminal))
# Plot the relationships among the eigengenes and the trait
sizeGrWindow(5,7.5);
par(cex = 0.9)
png("step7-Eigengene-dendrogram.png",width = 800,height = 600)
plotEigengeneNetworks(MET, "", marDendro = c(0,4,1,2), marHeatmap = c(3,4,1,2), cex.lab = 0.8, xLabelsAngle
= 90)
dev.off()
# Plot the dendrogram
sizeGrWindow(6,6);
par(cex = 1.0)
## 模塊的進(jìn)化樹
png("step7-Eigengene-dendrogram-hclust.png",width = 800,height = 600)
plotEigengeneNetworks(MET, "Eigengene dendrogram", marDendro = c(0,4,2,0),
plotHeatmaps = FALSE)
dev.off()
# Plot the heatmap matrix (note: this plot will overwrite the dendrogram plot)
par(cex = 1.0)
## 性狀與模塊熱
png("step7-Eigengene-adjacency-heatmap.png",width = 800,height = 600)
plotEigengeneNetworks(MET, "Eigengene adjacency heatmap", marHeatmap = c(3,4,2,2),
plotDendrograms = FALSE, xLabelsAngle = 90)
dev.off()
}
這個(gè)非常消耗計(jì)算資源和時(shí)間，所以建議選取其中部分基因作圖即可，我就沒有畫，而且根據(jù)下面的代碼選取部分基因來作圖！
然后隨機(jī)選取部分基因作圖
nSelect = 400
# For reproducibility, we set the random seed
set.seed(10);
select = sample(nGenes, size = nSelect);
selectTOM = dissTOM[select, select];
# There’s no simple way of restricting a clustering tree to a subset of genes, so we must re-cluster.
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;
TOMplot(plotDiss, selectTree, selectColors, main = "Network heatmap plot, selected genes")
這個(gè)圖湊數(shù)的意義居多，如果是把全部基因畫上去，可以很清楚的看到各個(gè)區(qū)塊顏色差異。
最后畫模塊和性狀的關(guān)系
# Recalculate module eigengenes
MEs = moduleEigengenes(datExpr, moduleColors)$eigengenes
## 只有連續(xù)型性狀才能只有計(jì)算
## 這里把是否屬于 Luminal 表型這個(gè)變量用0,1進(jìn)行數(shù)值化。
Luminal = as.data.frame(design[,3]);
names(Luminal) = "Luminal"
# Add the weight to existing module eigengenes
MET = orderMEs(cbind(MEs, Luminal))
# Plot the relationships among the eigengenes and the trait
sizeGrWindow(5,7.5);
par(cex = 0.9)
plotEigengeneNetworks(MET, "", marDendro = c(0,4,1,2), marHeatmap = c(3,4,1,2), cex.lab = 0.8, xLabelsAngle
= 90)
# 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, xLabelsAngle = 90)
step8:提取指定模塊的基因名
## step 8
# 主要是關(guān)心具體某個(gè)模塊內(nèi)部的基因
if(T){
# Select module
module = "brown";
# Select module probes
probes = colnames(datExpr) ## 我們例子里面的probe就是基因
inModule = (moduleColors==module);
modProbes = probes[inModule];
head(modProbes)
# 如果使用WGCNA包自帶的熱圖就很丑。
which.module="brown";
dat=datExpr[,moduleColors==which.module ]
plotMat(t(scale(dat)),nrgcols=30,rlabels=T,
clabels=T,rcols=which.module,
title=which.module )
datExpr[1:4,1:4]
dat=t(datExpr[,moduleColors==which.module ] )
library(pheatmap)
pheatmap(dat ,show_colnames =F,show_rownames = F) #對(duì)那些提取出來的1000個(gè)基因所在的每一行取出，組合起來為一個(gè)新的表達(dá)矩陣
n=t(scale(t(log(dat+1)))) # 'scale'可以對(duì)log-ratio數(shù)值進(jìn)行歸一化
n[n>2]=2
n[n< -2]= -2
n[1:4,1:4]
pheatmap(n,show_colnames =F,show_rownames = F)
group_list=datTraits$subtype
ac=data.frame(g=group_list)
rownames(ac)=colnames(n)
pheatmap(n,show_colnames =F,show_rownames = F,
annotation_col=ac )
# 可以很清晰的看到，所有的形狀相關(guān)的模塊基因
# 其實(shí)未必就不是差異表達(dá)基因。
}
有了基因信息，下游分析就很簡單了。包括GO/KEGG等功能數(shù)據(jù)庫的注釋
[圖片上傳失敗...(image-f522dd-1616404287721)]
Step9: 模塊的導(dǎo)出
主要模塊里面的基因直接的相互作用關(guān)系信息可以導(dǎo)出到cytoscape,VisANT等網(wǎng)絡(luò)可視化軟件。
# Recalculate topological overlap
TOM = TOMsimilarityFromExpr(datExpr, power = 6);
# Select module
module = "brown";
# Select module probes
probes = colnames(datExpr) ## 我們例子里面的probe就是基因名
inModule = (moduleColors==module);
modProbes = probes[inModule];
## 也是提取指定模塊的基因名
# Select the corresponding Topological Overlap
modTOM = TOM[inModule, inModule];
dimnames(modTOM) = list(modProbes, modProbes)
## 模塊對(duì)應(yīng)的基因關(guān)系矩陣
首先是導(dǎo)出到VisANT
vis = exportNetworkToVisANT(modTOM,
file = paste("VisANTInput-", module, ".txt", sep=""),
weighted = TRUE,
threshold = 0)
然后是導(dǎo)出到cytoscape
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.02,
nodeNames = modProbes,
nodeAttr = moduleColors[inModule]
);
如果模塊包含的基因太多，網(wǎng)絡(luò)太復(fù)雜，還可以進(jìn)行篩選，比如：
nTop = 30;
IMConn = softConnectivity(datExpr[, modProbes]);
top = (rank(-IMConn) <= nTop)
filter <- modTOM[top, top]
后面就是cytoscape自身的教程了，這里不再贅述，我博客有比較詳盡的介紹。
