3. 模块-临床信息关联并识别重要基因

3.1. 导入处理的数据：

首先，导入准备好的样本表达数据datExpr、临床信息datTraits以及构建的网络和模块networkConstruction-auto.RData：

# 设定工作目录
setwd("/home/cxt/data/R/")
# 导入包
library(WGCNA)
# 设置将strings不要转成factors
options(stringsAsFactors=F)

# 允许使用最大线程
allowWGCNAThreads()
# 或者直接指定线程数
# enableWGCNAThreads(nThreads=16)

# 导入表达与临床数据
lnames = load(file = "dataInput.RData");
lnames
# 导入共表达网络及模块
lnames = load(file = "networkConstruction-auto.RData");
lnames

3.2. 模块-临床信息关联定量

这步发分析主要是为例确定与我们感兴趣的临床特征有关的模块

在这项分析中，我们想要确定与定量与临床特征显著相关的模块。由于我们已经有了每个模块的信息(即，特征基因 Eigengene )，我们只是将特征基因与临床特征相关联，并寻找最重要的关联：

# 确定基因与样本个数
nGenes = ncol(datExpr);
nSamples = nrow(datExpr);
# 用颜色标签重新估计MEs（模块特征基因）
MEs0 = moduleEigengenes(datExpr, moduleColors)$eigengenes
MEs = orderMEs(MEs0)
moduleTraitCor = cor(MEs, datTraits, use = "p");
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples);

由于我们有相当多的模块和临床特征，适当的图形表示将有助于理解关联表。我们用相关值对每个关联进行颜色编码:

sizeGrWindow(10,6)
# 显示相关性及其p值
textMatrix =  paste(signif(moduleTraitCor, 2), "\n(",
                           signif(moduleTraitPvalue, 1), ")", sep = "");
dim(textMatrix) = dim(moduleTraitCor)
par(mar = c(6, 8.5, 3, 3));
# 在热图图中显示相关值
labeledHeatmap(Matrix = moduleTraitCor,
               xLabels = names(datTraits),
               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"))

从结果，我们可以确定了几个重要的模块-临床特征关联。这里我们把体重作为兴趣的特征

3.2. 基因与重要模块筛选（与临床相关）：Gene Significance and Module Membership

这步分析主要是为例计算每一个基因与临床特征的相关性，以及该基因与其所在模块的相关性

我们通过定义Gene Significance GS（绝对值）来量化个体基因与我们的兴趣特征(体重)之间的关联。

对于每个模块，我们还定义了一个Module Membership MM来量化模块特征基因（Eigengene ）和基因表达谱的相关性。

这使我们能够量化数组中所有基因与每个模块的相似性。

# 从临床信息中单独提取体重Weight的数据
weight = as.data.frame(datTraits$weight_g);
names(weight) = "weight"
# 获得模块的名称（颜色）
modNames = substring(names(MEs), 3)

## 计算Module Membership MM指标
geneModuleMembership = as.data.frame(cor(datExpr, MEs, use = "p"));
MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples));

names(geneModuleMembership) = paste("MM", modNames, sep="");
names(MMPvalue) = paste("p.MM", modNames, sep="");

## 计算Gene Significance GS指标
geneTraitSignificance = as.data.frame(cor(datExpr, weight, use = "p"));
GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples));

names(geneTraitSignificance) = paste("GS.", names(weight), sep="");
names(GSPvalue) = paste("p.GS.", names(weight), sep="");

3.3. 模块内分析：确定高MM和GS的基因

在这一步主要是为了展示模块中基因的GS和MM

基于GS和MM指标，我们可以筛选出一些重要的基因，这些基因即与体重高度相关，也在模块内很重要。比如，我们看一下与重量关联最大的棕色模块（MEbrown）。我们在MEbrown中绘制了基因重要性与模块中基因的关系的散点图：

## 获取棕色模块及其包含基因
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 body weight",
                   main = paste("Module membership vs. gene significance\n"),
                   cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)

如图所示，很明显，GS和MM是高度相关的，这说明与临床特征高度相关的基因通常也是与临床特征相关的模块中最重要的（中心）元素（基因）。

建议尝试使用其他重要临床特征/模块相关(例如，magenta, midnightblue, and red modules和体重)。

3.4. 保存网络分析的总结果

我们发现了与我们的兴趣特征高度相关的模块，并通过MM确定了模块内的核心元素（基因）。现在，我们将这些统计信息与基因注释合并在一起，并保存为一个文件，该文件总结了最重要的结果，可以在标准的电子表格软件(如MS)中进行检查：

我们的表达式数据只由探针ID名称注释

# 获得所有的探针ID（就是基因ID）
names(datExpr)

# 获得属于棕色模块的探针ID
names(datExpr)[moduleColors=="brown"]

为了便于解释结果，我们使用表达阵列制造商提供的探针注释文件将探针id连接到基因名称和公认的识别号码(Entrez代码)。

# 读取注释文件
annot = read.csv(file = "GeneAnnotation.csv");
dim(annot)
names(annot)

# 匹配表达矩阵中的探针文件
probes = names(datExpr)
probes2annot = match(probes, annot$substanceBXH)
# 检查是否有没有注释信息的探针
sum(is.na(probes2annot))
# 返回的结果应该是0

现在，我们为所有探针创建一个包含以下信息的Data Frame：probe ID、gene symbol、Locus Link ID (Entrez code)、module color、GS for weight 和 MM and p-values in all modules。模块将按照权重的重要性排序，最重要的模块在左边。

# 创建初始Data Frame
geneInfo0 = data.frame(substanceBXH = probes,
                      geneSymbol = annot$gene_symbol[probes2annot],
                      LocusLinkID = annot$LocusLinkID[probes2annot],
                      moduleColor = moduleColors,
                      geneTraitSignificance,
                      GSPvalue)
# 按与体重的相关性将模块重排
modOrder = order(-abs(cor(MEs, weight, use = "p")));
# 按选择的顺序添加MM信息
for (mod in 1:ncol(geneModuleMembership))
{
  oldNames = names(geneInfo0)
  geneInfo0 = data.frame(geneInfo0, geneModuleMembership[, modOrder[mod]], 
                         MMPvalue[, modOrder[mod]]);
  names(geneInfo0) = c(oldNames, paste("MM.", modNames[modOrder[mod]], sep=""),
                       paste("p.MM.", modNames[modOrder[mod]], sep=""))
}
# 首先按模块颜色对geneInfo变量中的基因排序，然后按geneTraitSignificance排序
geneOrder = order(geneInfo0$moduleColor, -abs(geneInfo0$GS.weight));
geneInfo = geneInfo0[geneOrder, ]

# 保存结果
write.csv(geneInfo, file = "geneInfo.csv")

3.5. 小总结

做完上述分析，在获得的结果中，我们就可以获得以下信息（比如我研究的是肿瘤）：

1. 与感兴趣临床特征（病人生存期）相关的模块
2. 获得模块中即与临床特征相关，也在模块中比较重要的基因
3. 对这些基因我们就可以进行下一步的分析