=====WGCNA实战（一）======
我们第一个实战采用的是官方提供的矩阵。
https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/Tutorials/
数据读入：
library(WGCNA)
library(reshape2)
library(stringr)
exprMat <- "LiverFemale3600.clean.txt"
trait <- " ClinicalTraits.csv"
options(stringsAsFactors = FALSE)
# 打开多线程
enableWGCNAThreads()
# 官方推荐"signed" 或 "signed hybrid"
type = "unsigned"
# 相关性计算
# 官方推荐 biweightmid-correlation & bicor
# corType: pearson or bicor
corType = "pearson"
corFnc =ifelse(corType=="pearson", cor, bicor)
# 对二元变量，如样本性状信息计算相关性时，
# 或基因表达严重依赖于疾病状态时，需设置下面参数
maxPOutliers =ifelse(corType=="pearson",1,0.05)
# 关联样品性状的二元变量时，设置
robustY =ifelse(corType=="pearson",T,F)
##导入数据##
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]
===数据筛选（可选）====
## 筛选中位绝对偏差前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)),]
## 转换为样品在行，基因在列的矩阵
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)，纵轴是无标度网络的评估参数，数值越高，网络越符合无标度特征(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")
===构建共表达网络====
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: 计算机能处理的最大模块的基因数量 (默认5000)；
# 4G内存电脑可处理8000-10000个，16G内存电脑可以处理2万个，32G内存电脑可以处理3万个， 计算资源允许的情况下最好放在一个block里面。
# corType: pearson or bicor
# numericLabels: 返回数字而不是颜色作为模块的名字，后面可以再转换为颜色
# saveTOMs：最耗费时间的计算，存储起来，供后续使用
# mergeCutHeight: 合并模块的阈值，越大模块越少；越小模块越多，冗余度越大；一般在0.15-0.3之间
# loadTOMs: 避免重复计算
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)
===共表达网络输出====
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矩阵，导出Cytoscape可用的数据===
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文件，可以直接导入cytoscape进行网络的可视化
//threshold 默认为0.5, 可以根据自己的需要调整，也可以都导出后在cytoscape中再调整
本文使用 文章同步助手 同步
网友评论