=====WGCNA實戰(zhàn)（一）======
我們第一個實戰(zhàn)采用的是官方提供的矩陣。
https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/Tutorials/
?
數(shù)據(jù)讀入：
library(WGCNA)
library(reshape2)
library(stringr)
exprMat <- "LiverFemale3600.clean.txt"??
trait <- " ClinicalTraits.csv"
options(stringsAsFactors = FALSE)
# 打開多線程
enableWGCNAThreads()
# 官方推薦"signed" 或 "signed hybrid"
type = "unsigned"
?
# 相關(guān)性計算
# 官方推薦 biweightmid-correlation & bicor
# corType: pearson or bicor
corType = "pearson"
?
corFnc =ifelse(corType=="pearson", cor, bicor)
# 對二元變量，如樣本性狀信息計算相關(guān)性時，
# 或基因表達嚴重依賴于疾病狀態(tài)時，需設(shè)置下面參數(shù)
maxPOutliers =ifelse(corType=="pearson",1,0.05)
?
# 關(guān)聯(lián)樣品性狀的二元變量時，設(shè)置
robustY =ifelse(corType=="pearson",T,F)
##導(dǎo)入數(shù)據(jù)##
dataExpr <- read.table(exprMat, sep='\t',row.names=1, header=T, quote="", comment="", check.names=F)
?
?
dim(dataExpr)
[1] 3600? 135
head(dataExpr)[,1:8]
===數(shù)據(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))
## 檢測缺失值
gsg = goodSamplesGenes(dataExpr, verbose =3)
?
if (!gsg$allOK)
{
? if(sum(!gsg$goodGenes)>0)
???printFlush(paste("Removing genes:",
????????????????????paste(names(dataExpr)[!gsg$goodGenes], collapse = ",")));
? if(sum(!gsg$goodSamples)>0)
??? printFlush(paste("Removingsamples:",
????????????????????paste(rownames(dataExpr)[!gsg$goodSamples], collapse = ",")));
? #Remove the offending genes and samples from the data:
?dataExpr = dataExpr[gsg$goodSamples, gsg$goodGenes]
}
===樣本層級聚類===
sampleTree = hclust(dist(dataExpr), method = "average")
plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="")
abline(h = 15, col = "red")
//從圖中可以看出有一個異常的sample F2_221，可以手動或者程序移除
clust = cutreeStatic(sampleTree, cutHeight = 15, minSize = 10)
table(clust)
keepSamples = (clust==1)
dataExpr2 = dataExpr[keepSamples, ]
dataExpr=dataExpr2
nGenes = ncol(dataExpr)
nSamples = nrow(dataExpr)
====確定軟閾值===
powers = c(c(1:10), seq(from = 12, to=20, by=2))
sft = pickSoftThreshold(dataExpr, powerVector = powers, verbose = 5)
sizeGrWindow(9, 5)
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")
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")
===構(gòu)建共表達網(wǎng)絡(luò)====
net = blockwiseModules(datExpr, power = 6,TOMType = "unsigned", minModuleSize = 30,reassignThreshold = 0, mergeCutHeight = 0.25,
numericLabels = TRUE, pamRespectsDendro = FALSE,saveTOMFileBase = "femaleMouseTOM",saveTOMs = TRUE,verbose = 3)
其中：
# power: 上一步計算的軟閾值?power = sft$powerEstimate
# maxBlockSize: 計算機能處理的最大模塊的基因數(shù)量 (默認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: 合并模塊的閾值，越大模塊越少；越小模塊越多，冗余度越大；一般在0.15-0.3之間
# loadTOMs: 避免重復(fù)計算
table(net$colors)
sizeGrWindow(12, 9)
mergedColors = labels2colors(net$colors)
plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],"Module colors",dendroLabels = FALSE, hang = 0.03,addGuide = TRUE, guideHang = 0.05)
moduleLabels = net$colors
moduleColors = labels2colors(moduleLabels)
dynamicColors <- labels2colors(net$unmergedColors)
plotDendroAndColors(net$dendrograms[[1]], cbind(dynamicColors,moduleColors),
? ? ? ? ? ? ? ? ? ? c("Dynamic Tree Cut", "Module colors"),
? ? ? ? ? ? ? ? ? ? dendroLabels = FALSE, hang = 0.5,
? ? ? ? ? ? ? ? ? ? addGuide = TRUE, guideHang = 0.05)
===共表達網(wǎng)絡(luò)輸出====
gene_module <-data.frame(ID=colnames(dataExpr), module=moduleColors)
gene_module =gene_module[order(gene_module$module),]
write.table(gene_module,file=paste0(exprMat,".gene_module.txt"),sep="\t",quote=F,row.names=F)? //我們也可以對每個module進行富集分析查看功能變化
//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)
MEs_colt = as.data.frame(t(MEs_col))
colnames(MEs_colt) = rownames(datExpr)
write.table(MEs_colt,file=paste0(exprMat,".module_eipgengene.V2.txt"),sep="\t",quote=F)
plotEigengeneNetworks(MEs_col, "Eigengene adjacency heatmap",?
? ? ? ? ? ? ? ? ? ? ? marDendro = c(3,3,2,4),
? ? ? ? ? ? ? ? ? ? ? marHeatmap = c(3,4,2,2), plotDendrograms = T,?
? ? ? ? ? ? ? ? ? ? ? xLabelsAngle = 90)
===篩選hub基因===
hubs = chooseTopHubInEachModule(dataExpr, colorh=moduleColors, power=power, type=type)
hubs
con <-nearestNeighborConnectivity(dataExpr, nNeighbors=50, power=power, type=type,corFnc = corType)
===獲取TOM矩陣，導(dǎo)出Cytoscape可用的數(shù)據(jù)===
load(net$TOMFiles[1], verbose=T)
TOM <- as.matrix(TOM)
?
dissTOM = 1-TOM
plotTOM = dissTOM^power
diag(plotTOM) = NA
probes = colnames(dataExpr)
dimnames(TOM) <- list(probes, probes)
cyt = exportNetworkToCytoscape(TOM,edgeFile = paste(exprMat, ".edges.txt", sep=""),nodeFile =paste(exprMat, ".nodes.txt", sep=""),weighted = TRUE,threshold = 0.01, nodeNames = probes, nodeAttr = moduleColors)??
//輸出node和edge文件，可以直接導(dǎo)入cytoscape進行網(wǎng)絡(luò)的可視化
//threshold 默認為0.5, 可以根據(jù)自己的需要調(diào)整，也可以都導(dǎo)出后在cytoscape中再調(diào)整
本文使用 文章同步助手 同步
