0.数据预整理

这是从GEO数据库下载的数据，GSE199335。

rm(list = ls())
library(tinyarray)
gse = "GSE199335"
geo = geo_download(gse,destdir=".",colon_remove = T)
geo$gpl

## [1] "GPL17400"

#View(geo$pd)
library(stringr)
Group = paste(geo$pd$genotype,geo$pd$age,sep="_") %>% 
  str_remove(" months of age| weeks of age") %>% 
  str_remove(" type") %>% 
  str_replace("/",".")
table(Group)

## Group
## R6.1_6 R6.2_9 wild_6 wild_9 
##      3      4      4      4

Group = factor(Group,levels = c("wild_6", "wild_9", "R6.1_6", "R6.2_9"))
ids <- AnnoProbe::idmap(geo$gpl,destdir = tempdir(),type = "soft")
ids = na.omit(ids)
exp = trans_array(geo$exp,ids,from = "ID")
pd = geo$pd
Group

##  [1] wild_6 wild_6 wild_6 R6.1_6 R6.1_6 wild_6 R6.1_6 wild_9 wild_9 wild_9
## [11] wild_9 R6.2_9 R6.2_9 R6.2_9 R6.2_9
## Levels: wild_6 wild_9 R6.1_6 R6.2_9

save(exp,Group,pd,file = "Dat.Rdata")

整理数据的代码，尤其是分组信息和探针注释，因数据而异。 下面的代码主要自WGCNA的官方文档，本强迫症进行了一些改动。

1.表达矩阵数据整理

首先需要有至少15个样本。 行是样本，列是基因。 不推荐使用全部基因，计算量太大 也不推荐使用差异分析所得的差异基因，包作者说的。 比较推荐的方法是按照方差或者mad去取前多少个基因，例如3500，5000，8000个，或者保留基因总数的前1/4。无标准答案。

rm(list = ls())
library(WGCNA)
library(tinyarray)
load("Dat.Rdata")
#exp = log2(geo$exp+1)
png("1.exp.png", width = 2000, height = 1200,res = 300)
boxplot(exp)
dev.off()

从这张图可以看出数据是取过log的，且没有异常样本，可以用。

datExpr0 = t(exp[order(apply(exp, 1, mad), decreasing = T)[1:5000],])
#datExpr0 = t(exp[order(apply(exp, 1, var), decreasing = T)[1:round(0.25*nrow(exp))],])

#datExpr0 = as.data.frame(t(exp))
datExpr0[1:4,1:4]

##              Bpifa6     Scd3     Myh4   Opn1sw
## GSM5970616 5.808989 5.964354 8.777054 7.488333
## GSM5970617 4.272349 4.216835 9.052218 7.843838
## GSM5970618 2.558643 4.095240 8.907886 7.464708
## GSM5970619 2.081443 4.061555 7.401055 2.912478

##              Bpifa6     Scd3     Myh4   Opn1sw
## GSM5970616 5.808989 5.964354 8.777054 7.488333
## GSM5970617 4.272349 4.216835 9.052218 7.843838
## GSM5970618 2.558643 4.095240 8.907886 7.464708
## GSM5970619 2.081443 4.061555 7.401055 2.912478
#rownames(datExpr0) = names(exp)[-c(1:8)]

1.1.基因过滤

主要看缺失值。GEO的芯片数据大多数没缺失值。

gsg = goodSamplesGenes(datExpr0, verbose = 3)

##  Flagging genes and samples with too many missing values...
##   ..step 1

##  Flagging genes and samples with too many missing values...
##   ..step 1
gsg$allOK # 返回TRUE则继续

## [1] TRUE

## [1] TRUE
if (!gsg$allOK){
  # 把含有缺失值的基因或样本打印出来
  if (sum(!gsg$goodGenes)>0)
    printFlush(paste("Removing genes:", paste(names(datExpr0)[!gsg$goodGenes], collapse = ", ")));
  if (sum(!gsg$goodSamples)>0)
    printFlush(paste("Removing samples:", paste(rownames(datExpr0)[!gsg$goodSamples], collapse = ", ")));
  # 去掉那些缺失值
  datExpr0 = datExpr0[gsg$goodSamples, gsg$goodGenes]
}

1.2.样本过滤

sampleTree = hclust(dist(datExpr0), method = "average")

png(file = "2.sampleClustering.png", width = 2000, height = 2000,res = 300)
par(cex = 0.6)
par(mar = c(0,4,2,0))
plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5,
     cex.axis = 1.5, cex.main = 2)
dev.off()

看有没有自己单独一个分支的样本， 如果有异常值就需要去掉，根据聚类图自己设置cutHeight 参数的值。

if(F){
  clust = cutreeStatic(sampleTree, cutHeight = 170, minSize = 10)
  table(clust) # 0代表切除的，1代表保留的
  keepSamples = (clust!=0)
  datExpr = datExpr0[keepSamples, ]
}
#没有异常样本就不需要去除
datExpr = datExpr0

2.表型信息整理

这个信息来自芯片数据的pData，也就是上面的pd。要求必须是数值型，要么像年龄那样的数字，要么搞成0，1，或者是1，2，3等。 这里利用了因子转数字的方法，转换成了数值。

library(stringr)
traitData = data.frame(genotype = as.numeric(as.factor(pd$genotype)),
                       age = as.numeric(as.factor(pd$age)))
str(traitData)

## 'data.frame':    15 obs. of  2 variables:
##  $ genotype: num  3 3 3 1 1 3 1 3 3 3 ...
##  $ age     : num  1 1 1 1 1 1 1 2 2 2 ...

## 'data.frame':    15 obs. of  2 variables:
##  $ genotype: num  3 3 3 1 1 3 1 3 3 3 ...
##  $ age     : num  1 1 1 1 1 1 1 2 2 2 ...
table(traitData[,1])

## 
## 1 2 3 
## 3 4 8

## 
## 1 2 3 
## 3 4 8
#pd
#dim(traitData)
names(traitData)
## [1] "genotype" "age"

这个数据有用的表型只有两列。

datTraits = traitData
sampleTree2 = hclust(dist(datExpr), method = "average")
# 用颜色表示表型在各个样本的表现: 白色表示低，红色为高，灰色为缺失
traitColors = numbers2colors(datTraits, signed = FALSE)
# 把样本聚类和表型绘制在一起
png(file = "3.Sample dendrogram and trait heatmap.png", width = 2000, height = 2000,res = 300)
plotDendroAndColors(sampleTree2, traitColors,
                    groupLabels = names(datTraits),
                    main = "Sample dendrogram and trait heatmap")
dev.off()

3.WGCNA正式开始

3.1 软阈值的筛选

设置一系列软阈值，范围是1-30之间，后面的数没必要全部画，就seq一下。

powers = c(1:10, seq(from = 12, to=30, by=2))
sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)

## pickSoftThreshold: will use block size 5000.
##  pickSoftThreshold: calculating connectivity for given powers...
##    ..working on genes 1 through 5000 of 5000
##    Power SFT.R.sq  slope truncated.R.sq mean.k. median.k. max.k.
## 1      1   0.4030  2.440          0.876 1590.00  1630.000 2170.0
## 2      2   0.0902  0.561          0.765  727.00   751.000 1230.0
## 3      3   0.0154 -0.169          0.673  396.00   402.000  785.0
## 4      4   0.1470 -0.444          0.720  240.00   236.000  538.0
## 5      5   0.3300 -0.651          0.777  156.00   147.000  394.0
## 6      6   0.4950 -0.813          0.850  107.00    96.200  304.0
## 7      7   0.6600 -0.916          0.952   76.90    65.000  245.0
## 8      8   0.7470 -1.100          0.980   56.90    45.100  209.0
## 9      9   0.8020 -1.270          0.993   43.30    32.200  186.0
## 10    10   0.8400 -1.420          0.992   33.70    23.600  168.0
## 11    12   0.8970 -1.600          0.975   21.60    13.400  141.0
## 12    14   0.9230 -1.660          0.954   14.70     8.010  121.0
## 13    16   0.9450 -1.670          0.952   10.40     5.000  106.0
## 14    18   0.9630 -1.660          0.959    7.70     3.210   94.5
## 15    20   0.9610 -1.630          0.951    5.86     2.130   84.7
## 16    22   0.9700 -1.590          0.961    4.58     1.450   76.5
## 17    24   0.9660 -1.560          0.958    3.65     1.020   69.5
## 18    26   0.9650 -1.520          0.958    2.97     0.724   63.5
## 19    28   0.9690 -1.490          0.966    2.45     0.513   58.2
## 20    30   0.9670 -1.450          0.967    2.05     0.378   53.6

sft$powerEstimate

## [1] 12

这个结果就是推荐的软阈值，拿到了可以直接用，无视下面的图。有的数据走到这一步会得到NA，也就是没得推荐。。。那就要看下面的图，选拐点。

根据我的经验，没有推荐软阈值，或者数字太大，后面跌跌撞撞走起来有些艰难哦，就得跑到前面重新调整表达矩阵里纳入的基因了。 cex1一般设置为0.9，不太合适（就是大多数软阈值对应的纵坐标都达不到0.9）时，可以设置为0.8或者0.85。一般不能再低了。

cex1 = 0.9
png(file = "4.Soft threshold.png", width = 2000, height = 1500,res = 300)
par(mfrow = c(1,2))
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=cex1,col="red")
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")
dev.off()

3.2 一步构建网络

如果前面没有得到推荐的软阈值，这里就要自己指定一个。根据上面的图，选择左图纵坐标第一个达到上面设置的cex1值的软阈值。

power = sft$powerEstimate
power

## [1] 12

net = blockwiseModules(datExpr, power = power,
                       TOMType = "unsigned", 
                       minModuleSize = 30, 
                       reassignThreshold = 0, 
                       mergeCutHeight = 0.25,
                       deepSplit = 2 ,
                       numericLabels = TRUE,
                       pamRespectsDendro = FALSE,
                       saveTOMs = TRUE,
                       saveTOMFileBase = "testTOM",
                       verbose = 3)

##  Calculating module eigengenes block-wise from all genes
##    Flagging genes and samples with too many missing values...
##     ..step 1
##  ..Working on block 1 .
##     TOM calculation: adjacency..
##     ..will not use multithreading.
##      Fraction of slow calculations: 0.000000
##     ..connectivity..
##     ..matrix multiplication (system BLAS)..
##     ..normalization..
##     ..done.
##    ..saving TOM for block 1 into file testTOM-block.1.RData
##  ....clustering..
##  ....detecting modules..
##  ....calculating module eigengenes..
##  ....checking kME in modules..
##      ..removing 1 genes from module 1 because their KME is too low.
##      ..removing 1 genes from module 2 because their KME is too low.
##      ..removing 1 genes from module 4 because their KME is too low.
##      ..removing 1 genes from module 6 because their KME is too low.
##  ..merging modules that are too close..
##      mergeCloseModules: Merging modules whose distance is less than 0.25
##        Calculating new MEs...

后面的结果如果不理想，比如划分的模块太少，或者青色、灰色的模块占到大多数，其他颜色都很少，或者颜色模块聚类是凌乱的，可以倒回来调整几个参数。

deepSplit 默认2，调整划分模块的敏感度，值越大，越敏感，得到的模块就越多； minModuleSize 默认30，参数设置最小模块的基因数，值越小，小的模块就会被保留下来； mergeCutHeight 默认0.25，设置合并相似性模块的距离，值越小，就越不容易被合并，保留下来的模块就越多。 https://zhuanlan.zhihu.com/p/34697561

可以试试调整，但是。。怎么说呢，主要看选择的基因是否给力，不给力的话调整了也就稍微好一点点。

不是很有必要去尝试分步法构建网络，得到的结果一样，可以调整的参数上面也都有。

此处展示得到了多少模块，每个模块里面有多少基因。

table(net$colors)

## 
##   0   1   2   3   4   5   6   7   8   9 
## 200 785 685 654 628 626 527 447 253 195

mergedColors = labels2colors(net$colors)
png(file = "5.DendroAndColors.png", width = 2000, height = 1200,res = 300)
plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
                    "Module colors",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)
dev.off()

好的结果就是每个颜色都差不多在一起，（青色和灰色不在考虑范围内）。然后青色和灰色的基因不要太多。

因为灰色代表没有合适的聚类，青色是基因数量的模块，比如你输入5000个基因，其中3000个都属于青色，剩下的模块基因数量太少，就很难受了。

3.3 保存每个模块对应的基因

moduleLabels = net$colors
moduleColors = labels2colors(net$colors)
MEs = net$MEs
geneTree = net$dendrograms[[1]]
gm = data.frame(net$colors)
gm$color = moduleColors
head(gm)

##        net.colors   color
## Bpifa6          7   black
## Scd3            7   black
## Myh4            9 magenta
## Opn1sw          8    pink
## Mb              9 magenta
## Myh13           9 magenta

genes = split(rownames(gm),gm$color)
save(genes,file = "genes.Rdata")

我这里把每个模块对应的基因存为了Rdata，用于数据挖掘下一步需求，提取基因。比如你需要某模块的基因与差异基因取交集等。

3.4 模块与表型的相关性

计算基因与表型的相关性矩阵

nGenes = ncol(datExpr)
nSamples = nrow(datExpr)
MEs0 = moduleEigengenes(datExpr, moduleColors)$eigengenes
MEs = orderMEs(MEs0)
moduleTraitCor = cor(MEs, datTraits, use = "p")
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)

#热图
png(file = "6.labeledHeatmap.png", width = 2000, height = 2000,res = 300)
# 设置热图上的文字（两行数字：第一行是模块与各种表型的相关系数；
# 第二行是p值）
# signif 取有效数字
textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
                   signif(moduleTraitPvalue, 1), ")", sep = "")
dim(textMatrix) = dim(moduleTraitCor)
par(mar = c(6, 8.5, 3, 3))
# 然后对moduleTraitCor画热图
labeledHeatmap(Matrix = moduleTraitCor,
               xLabels = names(datTraits),
               yLabels = names(MEs),
               ySymbols = names(MEs),
               colorLabels = FALSE,
               colors = blueWhiteRed (50),
               textMatrix = textMatrix,
               setStdMargins = FALSE,
               cex.text = 0.5,
               zlim = c(-1,1),
               main = paste("Module-trait relationships"))

dev.off()

我们希望找到和每个表型相关性较强的模块，正负相关都可。相关系数越大越好，如果能有个0.8，那结论就比较稳啦！没有的话0.6或0.7几也行，再小就不要拿来糊弄人了。

3.5. GS与MM

GS代表模块里的每个基因与形状的相关性

MM代表每个基因和所在模块之间的相关性，表示是否与模块的趋势一致。

modNames = substring(names(MEs), 3)
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="")

第几列的表型是最关心的，下面的i就设置为几。

与关心的表型相关性最高的模块赋值给下面的module。

i = 2
#module = "pink"
module = "turquoise"
assign(colnames(traitData)[i],traitData[i])
instrait = eval(parse(text = colnames(traitData)[i]))
geneTraitSignificance = as.data.frame(cor(datExpr, instrait, use = "p"))
GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples))
names(geneTraitSignificance) = paste("GS.", names(instrait), sep="")
names(GSPvalue) = paste("p.GS.", names(instrait), sep="")
png(file = paste0("7.MM-GS-scatterplot.png"), width = 2000, height = 2000,res = 300)
column = match(module, modNames) #找到目标模块所在列
moduleGenes = moduleColors==module #找到模块基因所在行
par(mfrow = c(1,1))
verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
                   abs(geneTraitSignificance[moduleGenes, 1]),
                   xlab = paste("Module Membership in", module, "module"),
                   ylab = "Gene significance",
                   main = paste("Module membership vs. gene significance\n"),
                   cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
dev.off()

我们希望看到强悍的正相关。数据挖掘里面，可以把整个模块的基因拿出来和别的步骤结果取交集了。

也可以通过这里找到GS和MM都大的基因。

f = data.frame(GS = abs(geneModuleMembership[moduleGenes, column]),
MM = abs(geneTraitSignificance[moduleGenes, 1]))
rownames(f) = rownames(gm[moduleGenes,])
head(f)

##                 GS        MM
## Rnu1b1   0.8667593 0.9482171
## ND6      0.8599513 0.9310548
## Vmn1r127 0.8839290 0.8336949
## Snora62  0.8862433 0.9287403
## Crygf    0.8759284 0.9884732
## Snord68  0.8733251 0.9027953

3.6.TOM

用基因相关性热图的方式展示加权网络，每行每列代表一个基因。 一般取400个基因画就够啦，拿全部基因去做电脑要烧起来了。

就想看看对角线附近红彤彤的小方块，以及同一个颜色基本都在一起，快乐。

nSelect = 400
set.seed(10)
dissTOM = 1-TOMsimilarityFromExpr(datExpr, power = 6)

## TOM calculation: adjacency..
## ..will not use multithreading.
##  Fraction of slow calculations: 0.000000
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.

select = sample(nGenes, size = nSelect)
selectTOM = dissTOM[select, select]
# 再计算基因之间的距离树(对于基因的子集，需要重新聚类)
selectTree = hclust(as.dist(selectTOM), method = "average")
selectColors = moduleColors[select]
library(gplots)
myheatcol = colorpanel(250,'red',"orange",'lemonchiffon')
png(file = "8.Sub400-netheatmap.png", width = 2000, height = 2000,res = 300)
plotDiss = selectTOM^7
diag(plotDiss) = NA #将对角线设成NA，在图形中显示为白色的点，更清晰显示趋势
TOMplot(plotDiss, selectTree, selectColors, col=myheatcol,main = "Network heatmap plot, selected genes")
dev.off()

3.7 模块与表型的相关性

MEs = moduleEigengenes(datExpr, moduleColors)$eigengenes
MET = orderMEs(cbind(MEs, instrait))
png(file = "9.Eigengene-dengro-heatmap.png", width = 2000, height = 2000,res = 300)
par(cex = 0.9)
plotEigengeneNetworks(MET, "", marDendro = c(0,4,1,2), marHeatmap = c(4,4,1,2), cex.lab = 0.8, xLabelsAngle
                      = 90)
dev.off()

也可以把上面的图分开来画

png(file = "10.Eigengene-dendrogram.png", width = 2000, height = 2000,res = 300)
par(cex = 1.0)
plotEigengeneNetworks(MET, "Eigengene dendrogram", marDendro = c(0,4,2,0),
                      plotHeatmaps = FALSE)
dev.off()

png(file = "11.Eigengene-heatmap.png",  width = 2000, height = 2000,res = 300)
par(cex = 1.0)
plotEigengeneNetworks(MET, "Eigengene adjacency heatmap", marHeatmap = c(4,5,2,2),
                      plotDendrograms = FALSE, xLabelsAngle = 90)
dev.off()