折腾了许久才复现出的WGCNA图-GSE106292你值得拥有

前些天老大布置了WGCNA的作业：下载GSE106292 数据集的 Excel表格如何读入R里面，做出作者文章中那样的图，但是分组很复杂，需要去文章中找细节内容，老大jimmy还写了篇推文：没有生物学背景的数据分析很危险 ，反反复复重做了好几次，终于做出和文章中的WGCNA的图和富集分析的结果非常接近的答案，还挺开心的。

要复现的图下

image-20191106102541372 image-20191106102656043 image-20191106102739038

1.下载整理数据

rm(list = ls())
options(stringsAsFactors = F)
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
options(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/")
library(GEOquery)
library(WGCNA)
f='GSE106292_eSet.Rdata'
# https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE106292

# 这个包需要注意两个配置，一般来说自动化的配置是足够的。
#Setting options('download.file.method.GEOquery'='auto')
#Setting options('GEOquery.inmemory.gpl'=FALSE)
if(!file.exists(f)){
  gset <- getGEO('GSE106292', destdir=".",
                 AnnotGPL = F,     ## 注释文件
                 getGPL = F)       ## 平台文件
  save(gset,file=f)   ## 保存到本地
}
load('GSE106292_eSet.Rdata')  ## 载入数据
class(gset)  #查看数据类型
length(gset)  #
class(gset[[1]])
# 因为这个GEO数据集只有一个GPL平台，所以下载到的是一个含有一个元素的list
a=gset[[1]] 
dat=exprs(a) 
dat<-read.csv('GSE106292_Human_Matrix_final.csv',sep = ',',row.names = 1)#rna-seq的数据以 Excel表格形式上传了，因此上一步通过dat=exprs(a)并不能获得表达矩阵信息，因此先下载表达矩阵再读进R里。
dim(dat)
dat<-dat[-1,]
rn_dat<-rownames(dat)
dat<-apply(dat,2,as.numeric)
dat=log2(dat+1)
dim(dat)
rownames(dat)<-rn_dat
dat[1:4,1:4]
colnames(dat)
colnames(dat)<-gsub('\\.',' ',colnames(dat))
colnames(dat)
pd=pData(a) #通过查看说明书知道取对象a里的临床信息用pData

原文描述WGCNA的段落是：
Here we implemented RNA sequencing to generate cell type-specific transcriptomes for chondrocytes, osteoblasts, myoblasts, tenocytes and ligamentocytes at 17 weeks post-conception (WPC) of human development. We then employed Weighted Gene Co-expression Network Analysis (WGCNA) to define tissue-specific gene modules that represent each cell type.

pd=pd[pd$`developmental stage:ch1`=="17 wks",]

原文描述tissue和临床信息表型对应关系的文字如下
chondrocytes from the knee, myoblasts from the quadriceps, endosteal osteoblasts from the femur, and ligamentocytes and tenocytes from the anterior and posterior cruciate ligament and Achilles tendon, respectively.

chondrocytes<-pd[pd$`tissue:ch1`=='Knee',]
myoblasts<-pd[pd$`tissue:ch1`=='quadriceps muscle',]
osteoblasts<-pd[pd$`tissue:ch1`=='Femur',]
ligamentocytes<-pd[pd$`tissue:ch1`=='Anterior/Posterior Cruciate Ligament',]
tenocytes<-pd[pd$`tissue:ch1`=='Achilles Tendon',]

接下来提取出与分组信息相对应的表达矩阵

pd1<-rbind(ligamentocytes,tenocytes,chondrocytes,myoblasts,osteoblasts)
pd2<-data.frame(row.names = rownames(pd1),
                title=pd1$title,
                tissue= rep(c('Ligamentocyte','Tenocyte','Chondrocyte','Myoblast','Osteoblast'),c(3,3,7,4,3)))
group_list<-rep(c('Ligamentocyte','Tenocyte','Chondrocyte','Myoblast','Osteoblast'),c(3,3,7,4,3))
dat1<-dat[,match(as.character(pd$title),colnames(dat))]
colnames(dat1)<-rownames(pd2)
dim(dat1)
dat1[1:4,1:4]
save(dat1,pd2,group_list,rn_dat,file = 'step1-input.Rdata')
load('step1-input.Rdata')

2.进行数据检查

rm(list = ls())
load('step1-input.Rdata')
dat<-dat1
dat<-apply(dat, 2, as.numeric)
table(apply(dat, 1, function(x) sum(x>1)))
table(apply(dat, 1, function(x) sum(x>1)==0))
table(apply(dat, 1, function(x) sum(x>1) > 5))

########################################  PCA  ##########################################


if(T){
  dat[1:4,1:4]
  dat<-log2(dat+1)
  ## 下面是画PCA的必须操作，需要看说明书。
  dat=t(dat)#画PCA图时要求是行名时样本名，列名时探针名，因此此时需要转换
  dat=as.data.frame(dat)#将matrix转换为data.frame
  dat=cbind(dat,group_list) #cbind横向追加，即将分组信息追加到最后一列
  #dat<-as.data.frame(dat)
  library("FactoMineR")#画主成分分析图需要加载这两个包
  library("factoextra") 
  # The variable group_list (index = 54676) is removed
  # before PCA analysis
  dat.pca <- PCA(dat[,-ncol(dat)], graph = FALSE)#现在dat最后一列是group_list，需要重新赋值给一个dat.pca,这个矩阵是不含有分组信息的
  fviz_pca_ind(dat.pca,
               geom.ind = "point", # show points only (nbut not "text")
               col.ind = dat$group_list, # color by groups
               palette = c("#9370DB", "#FF82AB", "#87CEFF", "#2E8B57", "#0000FF"),
               addEllipses = TRUE, # Concentration ellipses
               legend.title = "Groups")
  
  ggsave('PCA.png')
  
}

得到的下面的PCA图可以看到各组间还是分的很开的。

<img src="https://tva1.sinaimg.cn/large/006y8mN6gy1g96038tnwpj30me0im77a.jpg" alt="image-20191121215206122" style="zoom:50%;" />

3.热图-复现图a和图b：Top5000热图+各组相关性热图

#######################################  heatmap  ##########################################
rm(list = ls()) 
options(stringsAsFactors = F)
load(file = 'step1-input.Rdata')
dat<-dat1
cg=names(tail(sort(apply(dat,1,function(x){sum(x)})),5000))#apply按行（'1'是按行取，'2'是按列取），对每一行进行取表达量的最大值，从小到大排序，取最大的5000个
library(pheatmap)
pheatmap(dat[cg,],show_colnames =F,show_rownames = F) #对那些提取出来的1000个基因所在的每一行取出，组合起来为一个新的表达矩阵
n=t(scale(t(dat[cg,]))) # 'scale'可以对log-ratio数值进行归一化
n[n>3]=3
n[n< -3]= -3
n[1:4,1:4]
pheatmap(n,show_colnames =F,show_rownames = F)
ac=data.frame(g=group_list)
rownames(ac)=colnames(n) #把ac的行名给到n的列名，即对每一个探针标记上分组信息（是'noTNBC'还是'TNBC'）
pheatmap(n,show_colnames =T,show_rownames = F,annotation_names_col = F,
         annotation_col=ac,filename = 'heatmap_top5000.png')

共画了3张热图，最后一张热图展示如下图，与原文对比'Ligamentocyte'和'Chondrocyte'相比较其他组是高表达的。

image-20191121215218982

ac=data.frame(tissue=pd2$tissue)
rownames(ac)=colnames(dat)
M=cor(log(dat+1)) #计算列与例之间的相关性系数
pheatmap::pheatmap(M,annotation_row = ac,
                   annotation_col = ac,filename = 'cor.png') #这个地方很神奇， annotation_col是对列进行注释，那么出图结果是不仅加了列名，也加了行名

与原图对比，可以看到Ligamentocyte与Chondrocyte的相关性还是比较高的，同时Ligamentocyte与Tenocyte组的相关性也是相比较其他组高。而Myoblast组与其他组的相关性是最低的，与原文章中的也是相似的。当然，组内的各样本间相似性是最高的。

image-20191121215231252

4.WGCNA-复现图c：基因模块与性状相关性热图

参考：https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/faq.html

文章中关于进行WGCNA之前的数据为TPM数据(如下图所示)，根据一文看懂WGCNA 分析(2019更新版)：如果是芯片数据，那么常规的归一化矩阵即可，如果是转录组数据，最好是RPKM/TPM值或者其它归一化好的表达量。肉眼查看表达矩阵的数值大小，有的成百上千，有的为个位数甚至0，那么就需要用log2来进行归一化处理。所以直接用前面处理好的'step1-input.Rdata'就可以了，因为里面的dat1是已经经过log处理了的。

image-20191121215243030

接下来思考，我到底要不要过滤探针？
探针或基因可以通过平均表达或方差（或其鲁棒性强的MAD（中位数和中位数绝对偏差）进行过滤，因为低表达或不变好基因通常代表噪声。是否最好按平均表达式或方差进行筛选，这是一个争论的问题。两者都有优点和缺点，但更重要的是，它们倾向于筛选出相似的基因集，因为平均值和方差通常是相关的。

这篇文献复现时并没有采用通过基因表达量上的差异来过滤基因。经过查阅资料搜多到相关解释：WGCNA 被设计成一种无监督的分析方法，根据基因的表达特征对基因进行分组，通过基因表达量上的差异过滤后的基因，很可能就会导致形成一组相关基因就形成单个（或几个高度相关的）模块。我的理解是：如果通过基因表达量上的差异来过滤基因，就相当于类似人为地去划分模块了，而我们要的利用未经差异筛选后的表达矩阵来通过表达量高低与否将基因分在不同模块。

接下来再思考，那么如果不进行基因表达量上的差异来筛选基因，那么现在有20000多个基因，而且又有那么多表达量在很多样本中都为零的基因，我该如何过滤呢？见下图。

image-20191107014251551

经过搜索和尝试，我决定并不过滤很多基因，最后dat1<-dat1[!apply(dat1,1,function(x){sum(floor(x)==0)>15}),]来过滤基因，就是如果一个基因在20个样本中如果有超过15个样本表达量为0，那么这个基因我就不要了。过滤条件还是蛮低的，即使这么低，也还是过滤掉了20226-13178=7048个基因。

5.过滤基因，得到WGCNA的输入数据

rm(list = ls())
load('step1-input.Rdata')
dat1<-as.data.frame(dat1)
dat1[1:4,1:4]
apply(dat1,1,function(x){sum(floor(x)==0)>15})
dat1<-dat1[!apply(dat1,1,function(x){sum(floor(x)==0)>15}),]
dim(dat1)
apply(dat1,2,function(x){range(x)})
dat1<-t(dat1)
sampleTree = hclust(dist(dat1), method = "average")
png("step1_sampleClustering.png",width = 800,height = 600)
plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5,
     cex.axis = 1.5, cex.main = 2)
dev.off()
#得到接下来要进行WGCNA的输入数据
datExpr<-dat1
datTraits = data.frame(group_list=group_list) 
save(datExpr,datTraits,file = 'wgcna_input.Rdata')
load('wgcna_input.Rdata')

得到的样本聚类树可以看到，没有明显的离群样本，因此不需要剔除离群样本。

image-20191121215312541

6.获得Power值

rm(list = ls())
library(WGCNA)
load(file = 'wgcna_input.Rdata')
datExpr[1:4,1:4]

if(T){
  powers = c(c(1:10), seq(from = 12, to=30, by=2))
  # Call the network topology analysis function
  sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)
  #设置网络构建参数选择范围，计算无尺度分布拓扑矩阵
  png("step2-beta-value.png",width = 800,height = 600)
  # Plot the results:
  ##sizeGrWindow(9, 5)
  par(mfrow = c(1,2));
  cex1 = 0.7;
  # Scale-free topology fit index as a function of the soft-thresholding power
  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");
  # this line corresponds to using an R^2 cut-off of h
  abline(h=0.7,col="red")
  # Mean connectivity as a function of the soft-thresholding power
  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()
}
sft 
sft$powerEstimate
save(sft,file = "step2_beta_value.Rdata")

image-20191121215324327

查看软阈值，发现$powerEstimate这个函数计算出的推荐的软阈值是NA，但是通过搜索网上，发现这个问题同样也被问过，软阈值是可以自己挑选的，后面的步骤用自己挑选的就可以了。那么如何挑选？挑选SFT.R.sq的值尽量高，同时最大连通性mean.k.又不能太低。同时要根据下一步net$color生成的模块数目，我这里选择的power值是9,也就是R^2值为0.7。

image-20191121215333580

7.构建共表达矩阵

if(T){
  net = blockwiseModules(
    datExpr,
    power = 9,
    maxBlockSize = 6000,
    TOMType = "unsigned", minModuleSize = 30,
    reassignThreshold = 0, mergeCutHeight = 0.25,
    numericLabels = TRUE, pamRespectsDendro = FALSE,
    saveTOMs = F, 
    verbose = 3
  )
  table(net$colors) 
}

通过设置power=9，同时minModuleSize = 30(每个模块可包含的基因数目不能少于30个)，由于我看到有一些模块如从11到19，所包含的基因数目太少了，都低于100，所以我想在后面的代码中将minModuleSize设置为100。

image-20191107005850637

确定了贝塔值之后，就可以将关系矩阵转化为邻近矩阵，接下来就可以转换为tom重叠矩阵。为什么要转换为tom矩阵？是因为在wgcna中，认为模块是tom重叠性基因高的基因，所以需要计算基因和基因之间的tom重叠性，从而判断哪些基因应该属于同一个模块，哪些基因不在同一个模块。

8.模块可视化

if(T){
  #(1)网络构建 
  adjacency = adjacency(datExpr, power =9) #用softpower去计算邻接矩阵
  #(2) 邻近矩阵到拓扑矩阵的转换，Turn adjacency into topological overlap
  TOM = TOMsimilarity(adjacency);#邻接矩阵转换为tom重叠矩阵
  dissTOM = 1-TOM
  # (3) 聚类拓扑矩阵 Call the hierarchical clustering function
  geneTree = hclust(as.dist(dissTOM), method = "average");#将计算结果赋值给聚类树，计算方法是average，才是树枝
  # Plot the resulting clustering tree (dendrogram)
  sizeGrWindow(12,9)
  ## 这个时候的geneTree与一步法的 net$dendrograms[[1]] 性质类似，但是还需要进行进一步处理
  plot(geneTree, xlab="", sub="", main = "Gene clustering on TOM-based dissimilarity",
       labels = FALSE, hang = 0.04);
  #(4) 聚类分支的修整 dynamicTreeCut 
  # We like large modules, so we set the minimum module size relatively high:
  minModuleSize = 100;
  # Module identification using dynamic tree cut:
  dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM,
                              deepSplit = 2, pamRespectsDendro = FALSE,
                              minClusterSize = minModuleSize);
  table(dynamicMods)
  
  #绘画结果展示
  # Convert numeric lables into colors
  dynamicColors = labels2colors(dynamicMods)
  table(dynamicColors)
  # Plot the dendrogram and colors underneath
  #sizeGrWindow(8,6)
  png("dynamic tree cut.png",width = 800,height = 600)
  plotDendroAndColors(geneTree, dynamicColors, "Dynamic Tree Cut",
                      dendroLabels = FALSE, hang = 0.03,
                      addGuide = TRUE, guideHang = 0.05,
                      main = "Gene dendrogram and module colors")
  dev.off()
}

image-20191121215346305

##聚类结果相似模块的融合，Merging of modules whose expression profiles are very similar
#在聚类树中每一leaf是一个短线，代表一个基因，
#不同分之间靠的越近表示有高的共表达基因，将共表达极其相似的modules进行融合
# Calculate eigengenes
if(T){
  MEList = moduleEigengenes(datExpr, colors = dynamicColors)
  MEs = MEList$eigengenes
  # Calculate dissimilarity of module eigengenes
  MEDiss = 1-cor(MEs);#计算模块和模块之间的相关性和相异性
  # Cluster module eigengenes
  METree = hclust(as.dist(MEDiss), method = "average");
  # Plot the result
  #sizeGrWindow(7, 6)
  plot(METree, main = "Clustering of module eigengenes",
       xlab = "", sub = "")
  
  
  #建立基因模块后，可以将模块用颜色来区分，有些模块相似性高，就需要将模块合并。将模块特征基因进行聚类，在完成聚类后合并，0.15高度对应的相似度阈值就是0.85。具体的相似性阈值可以自行设置，进行聚类剪切后，就可以区分哪些模块相似性高，哪些模块相似性低，如下图。
  
  #选择有95%相关性的进行融合
  MEDissThres = 0.15#0.15剪切高度可修改  ####可以完成相似模块的合并，剪切高度是0.15，也就是将相似性高于0.85的模块进行了合并
  # Plot the cut line into the dendrogram
  abline(h=MEDissThres, col = "red")
  
  # Call an automatic merging function
  merge = mergeCloseModules(datExpr, dynamicColors, cutHeight = MEDissThres, verbose = 3)
  # The merged module colors
  mergedColors = merge$colors;
  # Eigengenes of the new merged modules:
  mergedMEs = merge$newMEs
  
  png("dynamicColors_mergedColors.png",width = 800,height = 600)
  plotDendroAndColors(geneTree, cbind(dynamicColors, mergedColors),
                      c("Dynamic Tree Cut", "Merged dynamic"),
                      dendroLabels = FALSE, hang = 0.03,
                      addGuide = TRUE, guideHang = 0.05)
  dev.off()
  
}

观察下图中哪些模块可以合并，设置融合线的高度。此处将融合高度设置为了0.15，完成相似模块的合并。剪切高度根据实际情况可修改。当剪切高度是0.15，也就是将相似性高于0.85的模块进行了合并。

image-20191121215355730 image-20191106212915144

经过融合后的基因模块与聚类树一同显示如下

image-20191121215412577

9.样本聚类可视化

if(T){
  nGenes = ncol(datExpr)
  nSamples = nrow(datExpr)
  #首先针对样本做个系统聚类
  datExpr_tree<-hclust(dist(datExpr), method = "average")
  #针对前面构造的样品矩阵添加对应颜色
  sample_colors1 <- numbers2colors(as.numeric(factor(datTraits$group_list)), 
                                   colors = c("green","blue","red","yellow","black"),signed = FALSE)
  
  ssss=as.matrix(data.frame(group_list=sample_colors1))                         
  par(mar = c(1,4,3,1),cex=0.8)
  png("sample-subtype-cluster.png",width = 800,height = 600)
  plotDendroAndColors(datExpr_tree, ssss,
                      groupLabels = colnames(sample),
                      cex.dendroLabels = 0.8,
                      marAll = c(1, 4, 3, 1),
                      cex.rowText = 0.01,
                      main = "Sample dendrogram and trait heatmap")
  dev.off()
}

image-20191121215422469

10.模块和性状的关系

## 这一步主要是针对于连续变量，如果是分类变量，需要转换成连续变量方可使用
table(datTraits)  
if(T){
  nGenes = ncol(datExpr)
  nSamples = nrow(datExpr)
  design1=model.matrix(~0+as.factor(datTraits$group_list))
  design=design1
  colnames(design) 
  colnames(design)=c(levels(as.factor(datTraits$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)
  
  sizeGrWindow(10,6)
  # Will display correlations and their p-values
  textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
                     signif(moduleTraitPvalue, 1), ")", sep = "");
  dim(textMatrix) = dim(moduleTraitCor)
  png("step5-Module-trait-relationships.png",width = 800,height = 1200,res = 120)
  par(mar = c(6, 8.5, 3, 3));
  # Display the correlation values within a heatmap plot
  labeledHeatmap(Matrix = moduleTraitCor,
                 xLabels = colnames(design),
                 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"))
  dev.off()
  table( labels2colors(net$colors))
}

<img src="https://tva1.sinaimg.cn/large/006y8mN6gy1g9605ywxb4j30nw0ygamj.jpg" alt="image-20191121215435698" style="zoom: 67%;" />

11.GO功能数据库的注释

table(moduleColors)
group_g=data.frame(gene=colnames(datExpr),
                   group=moduleColors)
save(group_g,file='wgcna_group_g.Rdata')



rm(list = ls())
load(file='wgcna_group_g.Rdata')
library(clusterProfiler)
# Convert gene ID into entrez genes
head(group_g)
tmp <- bitr(group_g$gene, fromType="SYMBOL", 
            toType="ENTREZID", 
            OrgDb="org.Hs.eg.db")

de_gene_clusters=merge(tmp,group_g,by.x='SYMBOL',by.y='gene')
table(de_gene_clusters$group)
head(de_gene_clusters)

list_de_gene_clusters <- split(de_gene_clusters$ENTREZID, 
                               de_gene_clusters$group)

library(ggplot2)
gcSample= list_de_gene_clusters
source('function.R')
# 下一步非常耗时，保守估计半个小时
# 主要是对我们的模块进行功能注释，就是GO/KEGG数据库
com_kegg_go(gcSample,'top5000')

• 对于Myoblasts在GO数据库中有富集通路如下

image-20191107022835232

• 对于Osteoblasts在GO数据库中有富集通路如下

image-20191107024052175

• 对于Chondrocyte在GO数据库中有富集通路如下

image-20191107024413247

上面的数据在得到的csv表格中同样能得到更多信息的验证。

最后友情宣传生信技能树

• 生物信息学“义诊”

• 生物信息学"拍卖会"

• 全国巡讲:R基础，Linux基础和RNA-seq实战演练 : 预告：12月28-30长沙站

• 广州珠江新城GEO数据挖掘滚动开班

• DNA及RNA甲基化数据分析与课题设计