学徒数据挖掘之IgG4-RD患者的WGCNA分析

学徒和学员已经陆续出师，是时候把生信技能树的舞台交给后辈了！

下面是《生信入门第7期》学员的分享

看完了⽣信技能树的系列WGCNA教程，见：

终于要开始自己的第⼀次实战笔记分享啦， 这次jimmy⽼师给了我⼀篇2020年⽂章的图表复现任务。

1.原文信息

原文标题：Multiple Processes May Involve in the IgG4-RD Pathogenesis: An Integrative Study via Proteomic and Transcriptomic Analysis

DOI： https://doi.org/10.3389/fimmu.2020.01795

分析思路：作者⽤了GSE40568和GSE66465两个数据集分别进行WGCNA分析，其中

• GSE40568为13个唇唾液腺组织样本，包括5个⼲燥综合症（SS）患者的样本、5个IgG4相关性疾病（IgG4-RD）患者的样本和3个正常（HC）样本；
• GSE66465为8个外周⾎单核细胞样本，包括两对治疗前后的gG4相关性疾病（IgG4-RD）患者的样本和4个对照的正常样本。

利⽤分组信息，通过WGCNA分析，筛选出和IgG4-RD最相关的⼀群基因，并进行KEGG分析以及下游验证。

WGCNA输入数据：变异系数>5%的基因

---

先看看原文的图~

Figure 1 是GSE40568数据集的WGCNA分析，图A是计算软阈值的结果图，这⾥作者选择的软阈值是20，图中红线的位置可能是0.75-0.8，图B是模块与性状关系的热图，C是GSE40568   中“绿松⽯”模块中的基因，D是模块和性状的关系图。Figure 2 则是GSE66465数据集的结果。

由于我不是免疫专业的，不清楚这个课题背景，因此不对分析结果做深⼊的解读。不过如果只看主干，略过其他地方，就能发现：作者的⽬的是通过WGCNA分析，找到与IgG4-RD最为相关的⼀群基因。

---

数据挖掘的核心是缩小目标基因

各种数据挖掘⽂章本质上都是要把目标基因集缩小，⽐如表达量矩阵通常是2万多个蛋⽩编码基因，不管是表达芯⽚还是RNA-seq测序的，采⽤何种程度的差异分析，最后都还有成百上千个⽬标基因。如果是临床队列，通常是会跟⽣存分析进⾏交集，或者多个数据集差异结果的交集，比如：多个数据集整合神器-RobustRankAggreg包，这样的基因集就是100个以内的数量了，但是仍然有缩小的空间，比如：lasso等统计学算法，最后搞成10个左右的基因组成signature即可顺利发表。

其实还有另外⼀个策略⽅向，有点类似于⼈⼯选择啦，通常是可以往热点靠，比如肿瘤免疫，相当于你不需要全部的两万多个基因的表达量矩阵进⾏后续分析，仅仅是拿着几千个免疫相关基因的表达矩阵即可。

最近比较热⻔的有：自噬基因，铁死亡，EMT基因，核受体基因家族，代谢基因。还有⼀个最搞笑的是m6a基因，完全是无厘头的基因集搞小，纯粹是为了搞小而搞小。

⽽我自己对GEO、TCGA数据库挖掘和WGCNA分析等的想法是这样的：（上面提到的几种）生信分析是不断划圈的过程，对一大群基因不断的划圈，缩小范围，最终得到我们想要的⼀⼩群基因（和分组标准最相关的⼀群基因）。

决定划圈（缩小范围）的因素，主要是分组（课题设计）和分析⽅法（例如WGCNA，差异分析，富集分析）。以下图为例：图中展示的是常见的差异分析和富集分析；图中的基因也可替换成蛋白、RNA等其他元素，划圈的方法也有很多，我就不⼀⼀举例了。

---

建⽴在划圈的理解上，这篇⽂章需要抓住的关键点就很清楚了：1.分组信息；2.分析方法；3.输⼊的数据和参数

• 分组标准，这⼀点可以在数据集里面找到相关信息；

• 分析方法—— WGCNA；

• 数据，作者⽤来⾃GSE40568和GSE66465两个数据集中变异系数（CV）> 5%的基因去做WGCNA分析；其中WGCNA分析的软阈值分别为20和12

搞清楚这几个关键点后，我们正式开始分析吧。

2.数据下载整理

# 这次我们还是使⽤GEOmirror下载数据，我只演示GSE40568数据集的代码和结果，
# GSE66465数据集是⼀样的过程
# remove global environment and all plots
rm(list = ls())
if (!is.null( dev.list() )) dev.off()
getwd()
# 实际上 jimmy⽼师很早以前就强调过，使⽤project形式组织⽂件，不要修改⼯作⽬录
setwd("D:/RCodingFile/biotree/GEO_WGCNA")
library(GEOquery)
library(GEOmirror)
eSet<- geoChina('GSE40568')
str(eSet)

# check the exprSet and creat group_list
eSet<- eSet[[1]]
exprSet <- exprs(eSet)
dim(exprSet)
head(exprSet[,1:4])
boxplot(exprSet,las=2)
exprSet

pdata<-pData(eSet)
table(pdata$source_name_ch1)
# 这⾥我点进去看了看pdata，前5个为SS样本，中间5个为IgG4RD样本，最后3个是正常样本
group_list = c(rep("SS",5), rep('IgG4RD',5),rep("HC",3))

# ⽤AnnoProbe注释探针
library(AnnoProbe)
gpl='GPL570'
probe2gene=idmap(gpl)
head(probe2gene)

# Delete unannotated probes
exprSet = exprSet[rownames(exprSet) %in% probe2gene[,1],]
probe2gene = probe2gene[match(rownames(exprSet), probe2gene[,1]),]
dim(exprSet)
dim(probe2gene)
tail(sort(table(probe2gene[,2])), n = 12L)

# Pick the most expressed gene
{
MAX = by(exprSet, probe2gene[,2],
function(x) rownames(x)[ which.max(rowMeans(x))])
MAX = as.character(MAX)
exprSet = exprSet[rownames(exprSet) %in% MAX,]
rownames( exprSet ) = probe2gene[ match( rownames( exprSet ), probe2gene[ , 1
] ), 2 ]
}
dim(exprSet)
exprSet[1:5,1:5]
AssayData = exprSet
save(AssayData,group_list, file = "./data/frist1step.Rdata")

基本上就是按照生信技能树的标准代码⼀步步下载数据。

3. 输入数据

# 上⼀步我们已经准备好了AssayData和group_list，这⼀步就是把数据输⼊到WGCNA
# 惯例清空环境
rm(list = ls())
if (!is.null( dev.list() )) dev.off()

load(file = "./data/frist1step.Rdata")
library("WGCNA")

# 根据分组情况修改列名
colnames( AssayData ) = paste( group_list, 1:ncol( AssayData ), sep = '_' )
dim( AssayData )
AssayData[1:5,]
# 求变异系数，我没找到合适的包就直接⾃⼰写了个函数
dat = AssayData
function2 <- function(dat){
    sd <- apply(dat,1,sd)
    mean <- apply(dat,1,mean)
    cv <- sd/mean
    return(cv[cv>=0.05])
}
cv = function2(dat)
AssayData = AssayData[names(cv),]
dim(AssayData)
# [1] 2270 13

咦，这⾥和原文不⼀致啊，原文是输⼊了变异系数＞5%的4364个基因，但是为什么我这里只有2270个基因的变异系数＞5% ？

经过不断摸索，我终于找到原因了~

eSet<- geoChina('GSE40568')
str(eSet)
eSet<- eSet[[1]]
exprSet <- exprs(eSet)
dim(exprSet)

dat = exprSet
function2 <- function(dat){
    sd <- apply(dat,1,sd)
    mean <- apply(dat,1,mean)
    cv <- sd/mean
    return(cv[cv>=0.05])
}

cv = function2(dat)
length(cv)
# [1] 4364

看，当不对原始数据进⾏任何处理时，CV＞5%的基因数刚好是4364。不同的⼈分析思路和代码习惯并不会完全相同，这里我选择剔除掉原始数据中没有注释基因的探针，并只保留相同基因中表达量最大的基因。（因为这些探针最后都没办法注释到基因，不可能去搞清楚到底它背后的⽣物学意义，还不如直接舍弃）

这样就可以避免剔除掉的无效数据（噪声）对后续WGCNA分析造成⼲扰，当然原文的思路也没有问题，只是每个人的习惯不⼀样而已。事实上，WGCNA的帮助⽂档也推荐剔除噪声，不过并没有指明需要剔除哪种噪声。

datExpr = t(AssayData)
dim( datExpr )
datExpr[1:4,1:4]
# ⾄此，WGCNA分析的数据都准备好了，开始绘制系统聚类树

4.系统聚类树

{
  nGenes = ncol(datExpr)
  nSamples = nrow(datExpr)
  datExpr_tree <- hclust(dist(datExpr), method = "average")
  par(mar = c(0,5,2,0))
  plot(datExpr_tree, main = "Sample clustering", sub = "", xlab = "",
       cex.axis = 0.9, cex.main = 1.2, cex.lab = 1, cex = 0.7)
  status_colors <- numbers2colors(as.numeric(factor(group_list)),
                                  colors = c("pink","darkgreen","lightblue"),signed = FALSE)
  par(mar = c(1,4,3,1),cex=0.8)
  plotDendroAndColors(datExpr_tree, status_colors,
                      groupLabels = "status",
                      cex.dendroLabels = 0.8,
                      marAll = c(1, 4, 3, 1),
                      cex.rowText = 0.01,
                      main = "Sample dendrogram and trait heatmap")
}

可以看到分组情况和聚类情况并不⼀致，不过这也很正常。临床组织样本中除了病变细胞，还有很多没有病变的间质细胞、血液淋巴和细胞外基质等，⼀致性不怎么好也很正常。

这也是经常细胞层面验证了的东西，换成临床样本就无法验证的原因。这个问题如果想解决，可以考虑取临床组织后消化为细胞，进行流式分选，挑出自己要的细胞再进行测序/芯片分析。不过这样得到的细胞量相对较少，操作也相对困难。

5.确定软阈值

powers = c(c(1:10), seq(from = 12, to = 20, by = 2))
## [1] 1 2 3 4 5 6 7 8 9 10 12 14 16 18 20
sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)

# Plot the results:
# Scale independence
par(mfrow = c(1,2))
cex1 = 0.75
plot(sft$fitIndices$Power,
     -sign(sft$fitIndices$slope) * sft$fitIndices$SFT.R.sq,
     xlab = "Soft Threshold (power)",
     ylab = "Scale Free Topology Model Fit,signed R^2",
     type = "n",
     main = paste("Scale independence"))
text(sft$fitIndices$Power,
     -sign(sft$fitIndices$slope) * sft$fitIndices$SFT.R.sq,
     labels = powers,
     cex = cex1,
     col = "red")
abline(h = 0.7, col = "RED")
# Mean connectivity
plot(sft$fitIndices$Power,
     sft$fitIndices$mean.k.,
     xlab = "Soft Threshold (power)",
     ylab = "Mean Connectivity",
     type = "n",
     main = paste("Mean connectivity"))
text(sft$fitIndices$Power,
     sft$fitIndices$mean.k.,
     labels = powers,
     cex = cex1,
     col = "red")
best_beta = sft$powerEstimate
best_beta
# 这里我选择的软阈值为20

6.构建共表达矩阵

net = blockwiseModules(datExpr, power = 20,
                       TOMType = "unsigned", minModuleSize = 30,
                       reassignThreshold = 0, mergeCutHeight = 0.15,
                       numericLabels = TRUE, pamRespectsDendro = FALSE,
                       saveTOMs = TRUE,
                       saveTOMFileBase = "femaleMouseTOM",
                       verbose = 3)


table(net$colors)
moduleColors <- labels2colors(net$colors)
# Recalculate MEs with color labels
MEs0 = moduleEigengenes(datExpr, moduleColors)$eigengenes
# 计算根据模块特征向量基因计算模块相异度：
MEDiss = 1 - cor(MEs0);
# Cluster module eigengenes
METree = hclust(as.dist(MEDiss), method = "average");
# Plot the result

plot(METree,
     main = "Clustering of module eigengenes",
     xlab = "",
     sub = "")
# 在聚类图中画出剪切线
MEDissThres = 0.25
abline(h = MEDissThres, col = "red")
merge_modules = mergeCloseModules(datExpr, moduleColors, cutHeight = MEDissThres, verbose
                                  = 3)
mergedColors = merge_modules$colors;
mergedMEs = merge_modules$newMEs;
plotDendroAndColors(net$dendrograms[[1]], cbind(moduleColors, mergedColors),
                    c("Dynamic Tree Cut", "Merged dynamic"),
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)

WGCNA把基因分为10个模块，其中灰色的模块较少，结果还是可以接受的。而原文中⼀共有13个模块，我的结果和原文不⼀致主要有两个原因：

• 输入的基因数量不同；

• 参数设定不同。

与模块大小相关的参数主要是blockwiseModules函数里面的minModuleSize、mergeCutHeight这两个参数。如果这两个参数越小，模块的大小也会越小，模块数量就会增多。我希望有较多的模块，所以设定为0.15，关于mergeCutHeight到底怎么设定并没有定论，我在网上找到⼀种说法大家可以参考参考

7.绘制TOM热图

# tom plot
nGenes = ncol(datExpr)
nSamples = nrow(datExpr)
geneTree = net$dendrograms[[1]]
dissTOM = 1 - TOMsimilarityFromExpr(datExpr, power = 20)
nSelect = 400
set.seed(10)
select = sample(nGenes, size = nSelect)
selectTOM = dissTOM[select, select]
selectTree = hclust(as.dist(selectTOM), method = "average")
selectColors = moduleColors[select]
sizeGrWindow(9,9)
plotDiss = selectTOM^7
diag(plotDiss) = NA
TOMplot(plotDiss, selectTree, selectColors, main = "Network heatmap plot, selected
genes")

8.模块与性状的关系

design = model.matrix(~0 + group_list)
colnames(design) = levels(group_list)
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)

# Will display correlations and their p-values
textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
                   signif(moduleTraitPvalue, 1), ")", sep = "");
dim(textMatrix) = dim(moduleTraitCor)
par(mar = c(6, 8.5, 3, 3));
# Display the correlation values within a heatmap plot
labeledHeatmap(Matrix = moduleTraitCor,
               xLabels = c("SS","IgG4RD","HC"),
               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"))

可以看到和IgG4-RD最相关的是turquoise模块，接下来会导出这个模块的基因进行富集分析。

MEs = moduleEigengenes(datExpr, moduleColors)$eigengenes
MET = orderMEs(MEs)
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,)

# names (colors) of the modules
modNames = substring(names(MEs), 3)
geneModuleMembership = as.data.frame(cor(datExpr, MEs, use = "p"));
## 算出每个模块跟基因的⽪尔森相关系数矩阵
## MEs是每个模块在每个样本⾥⾯的值
## datExpr是每个基因在每个样本的表达量
MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership),
nSamples));
names(geneModuleMembership) = paste("MM", modNames, sep="");
names(MMPvalue) = paste("p.MM", modNames, sep="");
geneTraitSignificance = as.data.frame(cor(datExpr,use = "p"));
GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance),
nSamples));
module = "turquoise"
column = match(module, modNames);
moduleGenes = moduleColors == module;
sizeGrWindow(7, 7);
par(mfrow = c(1,1));
verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),# MM > 0.8
                   abs(geneTraitSignificance[moduleGenes, 1]), # TS > 0.2
                   xlab = paste("Module Membership in", module, "module"),
                   ylab = "Gene significance for STATUS",
                   main = paste("Module membership vs. gene significance\n"),
                   cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)

9.挑选基因进行富集分析

# 原文中turquoise模块中共有934个基因，作者挑选其中edge-adjacency-value＞0.2的进行KEGG富集分析
# 我不是很清楚edge-adjacency-value到底指的是什么，所以我选择导出turquoise中 cor.geneTraitSignificance > 0.2的基因
# 我的处理方法得到的turquoise模块中有688个基因，挑出其中cor.geneTraitSignificance > 0.2的344个基因进行后续的富集分析
# kegg富集分析的代码请看以前的教程，这⾥我只放出结果

tmp1 = abs(geneTraitSignificance[moduleGenes, 1]) > 0.2
tmp2 = geneTraitSignificance[moduleGenes, ]
choose_gene = rownames(tmp2)[tmp1]

10.总结

辛苦各位小伙伴们看到这里了~

我只演示了GSE40568数据集的WGCNA分析结果，GSE66465的结果需要你自己动手去做哦。

由于不知道原作者的一些参数，我的结果会和原文不一致，并没有说谁对谁错~关键是要掌握WGCNA分析的目的和步骤。

WGCNA分析的目的：通过⼀定的算法，找到和性状（分组标准）最相关的一群基因。如果我们的分组标准是肿瘤与非肿瘤，突变与非突变，那我们得到的就是与肿瘤最相关的，或者是与突变最相关的⼀群基因。

WGCNA分析的步骤需要自己动手重复上面的代码才能掌握！注意某些函数参数的设定，遇到不清楚的可以查阅帮助文档或者谷歌寻找答案~

最后，我关于常见的生信数据挖掘的想法：依据分组标准，利用一定的计算方法（DEG、富集分析、WGCNA等），缩小范围，找到我们关心的⼀群（或者少数几个）基因~ 其中分组标准是核心，这个需要有⼀定的生物学背景；⽽计算方法则是重要的手段，需要提高自己的编程基础，两者兼顾才能做好自己的课题~

写在最后

本文是学徒数据挖掘专栏2020年度第一期，如果大家感兴趣2019年系列教程，欢迎查看：学徒数据挖掘专题半年目录汇总(生信菜鸟团周一见)