WGCNA仅仅是划分基因模块，其它都是附加分析

生信技能树

发布于 2023-09-04 15:04:19

1.4K00

代码可运行

文章被收录于专栏：生信技能树生信技能树

运行总次数：0

代码可运行

引言

曾老师给我分享了一篇数据挖掘的文章，里面的WGCNA非常奇怪，我之前没见过这样的模块与表型的相关性热图

看起来颜色、相关系数啥的非常对称，咋回事

看看实验设计：

GSE110993

Non-coding RNA profiling by high throughput sequencing

We used RNA sequencing to study expression changes of circulating miRNAs in a discovery sample of 20 IS patients and 20 matched healthy control subjects (HCs).

20个缺血性中风患者和20个正常对照外周血miRNA

但是这篇文章的作者使用样本的临床信息对实验进行了进一步分组

那就来做做看吧

获取数据

这篇文章属于数据挖掘，所以不是自己测的数据

查看使用的数据集，看起来也没直接提供表达矩阵

文章作者也是自己做的上游：

这一点需要注意，可以看到数据集GEO提供了上游处理流程和差异表达分析结果，我们拿这篇数据挖掘作者的过滤标准进行过滤得到的DEGs和数据集提供的并不一致（数量相差一半），而这篇文章中关于如何上游分析、差异表达分析、WGCNA的细节提的很少，但我们根据这一点可以初步判断：数据挖掘自己走的上游流程和数据集作者走的并不一致

另外作者这里行文前后的过滤标准不一致，也可能是我太菜了分辨不出来？我们这里使用后一种更宽松标准（|log2FoldChange|> 0.5, p-walue < 0.05）

后面经过师姐提醒，我在数据集发表文献的出版社网址上找到了miRNA的表达矩阵以及样本临床信息

RNA-Seq Identifies Circulating miR-125a-5p, miR-125b-5p, and miR-143-3p as Potential Biomarkers for Acute Ischemic Stroke

https://www.ahajournals.org/doi/suppl/10.1161/CIRCRESAHA.117.311572

这样我就不走上游了，miRNA测序定量的上游我还没走过，不知道和一般bulk一不一样，直接拿表达矩阵走WGCNA

带着问题出发

其实我在一开始接到曾老师给我的这个问题时，还不是非常熟悉WGCNA，只之前根据这两片推文跑了流程看看结果

这个WGCNA作业终于有学徒完成了！ RNA-seq入门实战（十一）：WGCNA加权基因共表达网络分析——关联基因模块与表型

在正式开始前我复习了这两个推文提到的问题以及流程代码，并检索学习了其他资料，有了一些自己的问题，贴在这里带着问题出发走后面的WGCNA，并希望能起到抛砖引玉的效果，得到小伙伴们在评论区的积极讨论

其它参考资料：

miRNA、LncRNA、CircRNA靠谱小结 搞懂RNA命名，miRNA、lncRNA、circRN不再傻傻分不清 转录组专题——WGCNA分析一问一答

问题：

Q1:作者使用的差异表达分析来进行的过滤，根据前面的学习并不推荐使用DEGs进行WGCNA，这会对结果产生什么样的影响？

RNA-seq入门实战（十一）：WGCNA加权基因共表达网络分析——关联基因模块与表型

Q2:两分组，小数量（几百）miRNA能做WGCNA吗？

转录组专题——WGCNA分析一问一答

相比于差异分析寻找在两组中显著差异表达的基因，趋势分析筛选与实验变化有关的最具代表性的基因集，WGCNA分析的独特优势体现在哪些呢？ 1）大样本数据挖掘：当有生物学重复时，样本数≥15个，不包含生物学重复时，样本数≥8个。 2）简洁归类：将大量的基因按照变化模式归类成不同的模块，简化整体分析难度。 3）注重基因关联：分析结果导入cytoscape中绘制基因共表达网络图，判断基因间潜在的调控关系，获得Hub基因。 4）关联表型维度：模块特征值与表型数据关联分析，关注与特定表型相关性高的模块。 ... WGCNA分析常见问题

小样本设计，如两组样本，每组三个生物学重复，能否开展WGCNA分析？

答：不能。主要有三个原因：（1）WGCNA分析是通过基因变化规律，将变化模式比较紧密的基因进行聚类生成不同的模块，表示它们之间存在潜在的调控关系。因此，按照基因变化模式进行聚类是WGCNA分析的关键。用于分析的样本需要包含非常丰富的变化信息，才能将基因归类成有生物学意义的模块。如果样本分组太少（如：两组、三组）变化模式比较弱，难以有效聚类；（2）WGCNA是以基因表达量相关系数为基础，得到基因变化模式的相关性。当样本数较少时，计算出的相关性无法判断哪些相关性是真实的，哪些相关性是随机波动产生的；（3）提出WCGNA分析的作者也推荐当样本是独立样本（不含生物学重复）时，建议样本数≥8，因为8种类型的样本，其变化模式已经十分丰富了。对于有生物学重复的样本，建议样本数≥15个。 ...

作者这里二分组还是拿来做了，感觉这里就没利用到WGCNA相较于差异表达分析用来处理多分组的优点，只利用到了基因按照变化模式归类成不同的模块

WGCNA一般流程

1.读取数据

数据集作者提供了三种表达矩阵

分别是未处理的counts矩阵、以及DESeq2均一化矩阵、EdgeR处理的TMM矩阵

#####读取数据#######
rm(list=ls())
library(tidyverse)
library(edgeR)
library(WGCNA)
library(FactoMineR)
library(factoextra)  
library(tidyverse) # ggplot2 stringer dplyr tidyr readr purrr  tibble forcats
library(data.table) #多核读取文件
# cpm表达矩阵
miRNA <- read.table("miRNA_counts.csv",header = TRUE,sep = ",")
miRNA <- column_to_rownames(miRNA,"X")
# cpm
range(miRNA)  # [1]       0 1920696
express_cpm <- log2(cpm(miRNA)+ 1)  # cpm {edgeR}
range(express_cpm)  # [1]  0.00000 17.75535
# 筛选MAD前1000的miRNA 不做差异分析
# 20000-5000 4000-1000
keep <- express_cpm[order(apply(express_cpm,1,mad), decreasing = T)[1:2000],]
keep_data <- express_cpm[rownames(keep),]
# 转置
dim(keep_data)  # [1] 1000   40
datExpr0 <- as.data.frame(t(keep_data))

# 分组
datTraits <- data.frame(
  row.names = rownames(datExpr0),
  exp_group = rep(c("IS","HC"),each=20),
  exp_group_No = rep(c(1,2),each=20)
)

2.PCA和聚类树

############PCA和聚类树图##########

### 绘制样品的系统聚类树
if(T){
  # 对处理完后的表达矩阵进行聚类，得到样本聚类树
  sampleTree <- hclust(dist(datExpr0), method = "average")
  par(mar = c(0,5,2,0))
  # 绘制样品的系统聚类树，可以查看样本是否聚类合理
  plot(sampleTree, main = "Sample clustering", sub="", xlab="", cex.lab = 2, 
       cex.axis = 1, cex.main = 1,cex.lab=1)
  ## 若样本有性状、表型，可以添加对应颜色，查看是否聚类合理
  # numbers2colors将样本的trait信息转化为对应的颜色，用于在样品聚类树上对照查看
  sample_colors <- numbers2colors(as.numeric(factor(datTraits$exp_group)), 
                                  colors = rainbow(length(table(datTraits$exp_group))), 
                                  signed = FALSE)
  ## 绘制样品的系统聚类树及对应性状
  par(mar = c(1,4,3,1),cex=0.8)
  pdf("step1_Sample dendrogram and trait.pdf",width = 8,height = 6)
  # 结合样品聚类树和性状信息，绘制样品的聚类树及对应性状
  plotDendroAndColors(sampleTree, sample_colors,
                      groupLabels = "trait",
                      cex.dendroLabels = 0.8,
                      marAll = c(1, 4, 3, 1),
                      cex.rowText = 0.01,
                      main = "Sample dendrogram and trait" )
  dev.off()
}

##若存在显著离群点；剔除掉（联系后面绘图解读）
if(F){
  # 检查样品的异常点，将离群值样品移除
  clust <- cutreeStatic(sampleTree, cutHeight = 23500, minSize = 10) # cutHeight根据实际情况而定
  table(clust)
  keepSamples <- (clust==1)
  datExpr0 <- datExpr0[keepSamples, ]
  datTraits <- datTraits[keepSamples,]
  dim(datExpr0) 
}

### 判断数据质量 : PCA进行分组查看
# 使用PCA对表达矩阵进行降维，便于后面的聚类和可视化
group_list <- datTraits$exp_group
dat.pca <- PCA(datExpr0, graph = F)  # 作者是datExpr我加一个0
# 绘制PCA的结果，用于查看数据是否存在分类趋势
pca <- fviz_pca_ind(dat.pca,
                    title = "Principal Component Analysis",
                    legend.title = "Groups",
                    geom.ind = c("point","text"), #"point","text"
                    pointsize = 2,
                    labelsize = 4,
                    repel = TRUE, #标签不重叠
                    col.ind = group_list, # 分组上色
                    axes.linetype=NA,  # remove axeslines
                    mean.point=F#去除分组中心点
) +
  theme(legend.position = "none")+  # "none" REMOVE legend
  coord_fixed(ratio = 1) #坐标轴的纵横比
pca
ggsave(pca,filename= "step1_Sample PCA analysis.pdf", width = 8, height = 8)

##
# PCA和聚类树图都可看出此次聚类效果不好的
# 并且不是靠剔除离群值就可以分开

# 迫不得已 跟作者一样用DEGs试试

我们在代码中对前面keep变量进行变换，可以发现MAD前1000和前2000的miRNA拿来聚类分组效果都不是很好

那就试试DEGs吧：

#######+DEGs###########

# 默认counts所以要求整数 内部有log2
#1 查看分组信息和表达矩阵数据
# keep <- rowSums(miRNA>0) >= floor(0.75*ncol(miRNA))
# filter_count <- miRNA[keep,] #获得filter_count矩阵
# exprSet <- filter_count
# dim(exprSet)
# ## keep <- rowSums(miRNA>0) >= floor(0.75*ncol(miRNA)) 
# ## 常规过滤去太多
# ## [1] 385  40
#
#不过滤看看 反正DESeq2有一套自过滤
exprSet <- miRNA
table(group_list)
## group_list
# HC IS 
# 20 20 
# 加载包
library(DESeq2)
#2 第一步，构建DESeq2的DESeq对象dds(建立DEseq数据矩阵)
colData <- data.frame(row.names=colnames(exprSet),group_list=group_list)

dds <- DESeqDataSetFromMatrix(countData = exprSet,
                              colData = colData,
                              design = ~ group_list)
# > dds
# class: DESeqDataSet 
# dim: 4446 40 
# metadata(1): version
# assays(1): counts
# rownames(4446): hsa-let-7a-1 hsa-let-7a-2 ... hsa-miR-99b-3p
# hsa-miR-99b-5p
# rowData names(0):
#   colnames(40): LH_41_hg19 LH_42_hg19 ... LH_79_hg19 LH_80_hg19
# colData names(1): group_list

#3 第二步，进行差异表达分析
dds2 <- DESeq(dds)
#4 提取差异分析结果，trt组对untrt组的差异分析结果
tmp <- results(dds2,contrast=c("group_list","IS","HC"))
DEG_DESeq2 <- as.data.frame(tmp[order(tmp$padj),])
head(DEG_DESeq2)  # 联系之前的内容DEseq2和edge会存在p值为0的情况
# baseMean log2FoldChange     lfcSE      stat       pvalue
# hsa-miR-1        14.25921      -2.507691 0.4733685 -5.297544 1.173706e-07
# hsa-mir-1-1      13.14767      -2.876105 0.5401860 -5.324287 1.013498e-07
# hsa-miR-3184-5p 201.71927      -1.241052 0.2651975 -4.679730 2.872532e-06
# hsa-miR-423-3p  201.71927      -1.241052 0.2651975 -4.679730 2.872532e-06
# hsa-mir-660      24.50594      -1.696137 0.3747271 -4.526327 6.001761e-06
# hsa-miR-660-5p   24.46895      -1.693303 0.3746032 -4.520258 6.176421e-06
# padj
# hsa-miR-1       2.253515e-05
# hsa-mir-1-1     2.253515e-05
# hsa-miR-3184-5p 2.757631e-04
# hsa-miR-423-3p  2.757631e-04
# hsa-mir-660     3.952910e-04
# hsa-miR-660-5p  3.952910e-04
# 去除差异分析结果中包含NA值的行
DEG_DESeq2 = na.omit(DEG_DESeq2)
dim(DEG_DESeq2)  # [1] 384   6

作者用的就是DEGs，这里看看我们前面提到了作者行文中过滤阈值不清楚的地方：

发现还是只能按照后一种较为宽松的标准来做

同时可以发现，数据集作者提供的差异表达基因数量和我们拿他提供的表达矩阵进行差异分析得到的差不多，说明这个表达矩阵差异分析过程还是可靠的

但100多个DEGs和这篇数据挖掘文章作者说的差别很大，应该是作者自己走了上游分析的原因

想看看文章作者到底用多少miRNA做WGCNA

但是这表格目前处于无效状态，网站里打不开

所以又陷入困境

使用MAD选取差异基因发现没办法很好地将实验分组分开

我猜文章作者也是这个原因所以拿DEGs做WGCNA

但是作者文章关于到底多少DEGs走了WGCNA，文章表达较为矛盾，也找不到过程文件

死马当做活马医看看这我们得到的100多个差异表达miRNA能不能很好区分实验分组，继续往下走，本文主要是看看作者那个模块与表型的相关性热图是怎么回事：

# 筛选上下调差异基因，设定阈值
fc_cutoff <-1
fdr <- 0.05#p值
DEG_DESeq2$regulated <- "normal"
loc_up <- intersect(which(DEG_DESeq2$log2FoldChange>log2(fc_cutoff)),
                    which(DEG_DESeq2$pvalue<fdr))
loc_down <- intersect(which(DEG_DESeq2$log2FoldChange< (-log2(fc_cutoff))),
                      which(DEG_DESeq2$pvalue<fdr))
DEG_DESeq2$regulated[loc_up] <- "up"
DEG_DESeq2$regulated[loc_down] <- "down"
table(DEG_DESeq2$regulated)
# 换阈值试试
# down normal     up 
# 110    262     12

# fc_cutoff <-2.0
# fdr <- 0.05#p值
# up :padj<fdr log2FoldChange>log2(fc_cutoff))
# down normal     up 
# 36    342      6 
# DEGs少得可怜

### 绘制样品的系统聚类树
DEGs <- DEG_DESeq2 %>% filter(regulated != "normal") %>% rownames()
DEGs_df <-as.data.frame(t(express_cpm[rownames(express_cpm) %in% DEGs,]))
if(T){
  # 对处理完后的表达矩阵进行聚类，得到样本聚类树
  sampleTree <- hclust(dist(DEGs_df), method = "average")
  par(mar = c(0,5,2,0))
  # 绘制样品的系统聚类树，可以查看样本是否聚类合理
  plot(sampleTree, main = "Sample clustering", sub="", xlab="", cex.lab = 2, 
       cex.axis = 1, cex.main = 1,cex.lab=1)
  ## 若样本有性状、表型，可以添加对应颜色，查看是否聚类合理
  # numbers2colors将样本的trait信息转化为对应的颜色，用于在样品聚类树上对照查看
  sample_colors <- numbers2colors(as.numeric(factor(datTraits$exp_group)), 
                                  colors = rainbow(length(table(datTraits$exp_group))), 
                                  signed = FALSE)
  ## 绘制样品的系统聚类树及对应性状
  par(mar = c(1,4,3,1),cex=0.8)
  pdf("step1_Sample dendrogram and trait.pdf",width = 8,height = 6)
  # 结合样品聚类树和性状信息，绘制样品的聚类树及对应性状
  plotDendroAndColors(sampleTree, sample_colors,
                      groupLabels = "trait",
                      cex.dendroLabels = 0.8,
                      marAll = c(1, 4, 3, 1),
                      cex.rowText = 0.01,
                      main = "Sample dendrogram and trait" )
  dev.off()
}

# 使用PCA对表达矩阵进行降维，便于后面的聚类和可视化
group_list <- datTraits$exp_group
dat.pca <- PCA(DEGs_df, graph = F)  # 作者是datExpr我加一个0
# 绘制PCA的结果，用于查看数据是否存在分类趋势
pca <- fviz_pca_ind(dat.pca,
                    title = "Principal Component Analysis",
                    legend.title = "Groups",
                    geom.ind = c("point","text"), #"point","text"
                    pointsize = 2,
                    labelsize = 4,
                    repel = TRUE, #标签不重叠
                    col.ind = group_list, # 分组上色
                    axes.linetype=NA,  # remove axeslines
                    mean.point=F#去除分组中心点
) +
  theme(legend.position = "none")+  # "none" REMOVE legend
  coord_fixed(ratio = 1) #坐标轴的纵横比
pca
ggsave(pca,filename= "step1_Sample PCA analysis.pdf", width = 8, height = 8)

##保存数据
datExpr <- DEGs_df
nGenes <- ncol(DEGs_df)
nSamples <- nrow(DEGs_df)
save(nGenes,nSamples,datExpr,datTraits,file="step1_input.Rdata")

看起来比之前效果好

3.挑选最佳阈值power

############################### 挑选最佳阈值power ###################################
# 挑选最佳阈值power
rm(list = ls())  
load("step1_input.Rdata")
# 设置一个R平方值的截断点，作为判断拟合优度的标准，默认值为0.8
R.sq_cutoff = 0.8  #设置R^2 cut-off
if(T){
  # Call the network topology analysis function
  #设置power参数选择范围 1到20，每次增加1，以及22到30，每次增加2
  # 为了在选择最佳power时，尽可能涵盖可能的取值范围
  powers <- c(seq(1,20,by = 1), seq(22,30,by = 2)) 
  # 用来选择最佳power参数的，
  # 输入参数包括表达数据datExpr、网络类型（这里是无向网络）、power取值的向量、R平方的截断值、以及是否输出过程信息
  # 输出结果为一个包含有关于不同power下拟合优度、均值连接度、平均k值等的信息表
  sft <- pickSoftThreshold(datExpr, 
                           networkType = "unsigned",
                           powerVector = powers, 
                           RsquaredCut = R.sq_cutoff,  
                           verbose = 5)
  #SFT.R.sq > 0.8 , slope ≈ -1
  # 将拟合优度图和均值连接度图保存至PDF文件中
  pdf("step2_power-value.pdf",width = 16,height = 12)
  # Plot the results: 寻找拐点，确认最终power取值
  par(mfrow = c(1,2));
  cex1 = 0.9;
  
  # Scale-free topology fit index as a function of the soft-thresholding power
  # 画出在不同power下拟合优度的趋势图，x轴为power，y轴为拟合优度（signed R2值）
  # 拟合优度要乘以正负1，这是因为负的拟合优度意味着模型拟合得更差
  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")
  # 在每个数据点上添加power的数值标注，颜色为红色
  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
  # 画一条水平的线，标注R平方的截断点，颜色为红色
  abline(h=R.sq_cutoff ,col="red")
  
  # Mean connectivity as a function of the soft-thresholding power
  # 画出在不同power下均值连接度的趋势图，x轴为power，y轴为均值连接度
  plot(sft$fitIndices[,1], sft$fitIndices[,5],
       xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n")
  # 在每个数据点上添加power的数值标注，颜色为红色
  text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")
  # 画一条水平的线，标注均值连接度的截断点，均值连接度大于100的边自动被过滤掉，颜色为红色
  abline(h=100,col="red")
  dev.off()
}  # 可看出拐点大致在power为16时出现，且各项数值基本满足挑选标准，因此设定power=16

# 打印出在不同power下估计的最佳power取值
sft$powerEstimate  #查看估计的最佳power [1] NA
power = sft$powerEstimate

# 若无向网络在power小于15或有向网络power小于30内，没有一个power值使
# 无标度网络图谱结构R^2达到0.8且平均连接度在100以下，可能是由于
# 部分样品与其他样品差别太大。这可能由批次效应、样品异质性或实验条件对
# 表达影响太大等造成。可以通过绘制样品聚类查看分组信息和有无异常样品。
# 如果这确实是由有意义的生物变化引起的，也可以使用下面的经验power值。

# 如果最佳power取值为空
# 根据数据样本数确定经验power取值：
# 如果样本数小于20，则unsigned类型网格选择power=9，否则选择power=等一系列取值；
# 如果样本数大于等于20，则unsigned类型网格选择power=6，否则选择power=12。这个默认给出较为通用的power取值方案，但是具体的数据情况下是否适用需再加以考虑。
if(is.na(power)){
  # 官方推荐 "signed" 或 "signed hybrid"
  # 为与原文档一致，故未修改
  type = "unsigned"
  nSamples=nrow(datExpr)
  power = ifelse(nSamples<20, ifelse(type == "unsigned", 9, 18),
                 ifelse(nSamples<30, ifelse(type == "unsigned", 8, 16),
                        ifelse(nSamples<40, ifelse(type == "unsigned", 7, 14),
                               ifelse(type == "unsigned", 6, 12))  
                 )
  )
}
# 将选择好的最佳power值和soft-thresholding的计算结果保存至RData文件中
save(sft, power, file = "step2_power_value.Rdata")

这里脚本选取的为经验power 6

4.一步法构建加权共表达网络，识别基因模块

######一步法构建加权共表达网络，识别基因模块######

# 主要调整参数为minModuleSize和mergeCutHeight , 数值越大模块越少
rm(list=ls())
load(file = "step1_input.Rdata")
load(file = "step2_power_value.Rdata")
if(T){
  cor <- WGCNA::cor
  # blockwiseModules函数构建加权共表达网络，并识别基因模块
  net <- blockwiseModules(
    datExpr,  # datExpr是包含基因表达数据的矩阵
    power = power,  # power是上一步计算得到的软阈值幂次
    maxBlockSize = ncol(datExpr),  # 计算机能处理的最大模块的基因数量
    corType = "pearson", #计算相关性的方法默认为"pearson","bicor"则更能考虑离群点的影响
    networkType = "unsigned",  # 无向网络 
    TOMType = "unsigned",   # 建立邻接矩阵和TOM矩阵时，是否考虑正负相关性，unsigned表示只考虑正相关性
    # minModuleSize指每个模块中基因的最少数目
    minModuleSize = 10,    ##越大模块越少
    # mergeCutHeight指合并模块的阈值
    mergeCutHeight = 0.25, ##越大模块越少
    numericLabels = TRUE,  # 返回数字作为模块的名称，后面可以再转换为颜色
    saveTOMs = F,  # 是否存储TOM矩阵，TOM矩阵计算最耗费时间的步骤之一，如果需要多次分析，可以将其存储起来供后续使用
    verbose = 3  # 控制输出信息的详细程度，数值越大输出的信息越多
  )
  # 统计每个模块中基因的数量
  table(net$colors) 
  # 0  1  2  3  4 
  # 18 56 22 16 10
  # power: 上一步计算的软阈值
  # maxBlockSize:计算机能处理的最大模块的基因数量(默认5000),16G内存可以处理2万个，
  # 计算资源允许的情况下最好放在一个block里面。
  # corType：计算相关性的方法；可选pearson(默认)，bicor。后者更能考虑离群点的影响。
  # networkType:计算邻接矩阵时，是否考虑正负相关性；默认为"unsigned",可选"signed", "signed hybrid"
  # TOMType：计算TOM矩阵时，是否考虑正负相关性；默认为"signed",可选"unsigned"。但是根据幂律转换的邻接矩阵(权重)的非负性，所以认为这里选择"signed"也没有太多的意义。
  # numericLabels: 返回数字而不是颜色作为模块的名字，后面可以再转换为颜色
  # saveTOMs：最耗费时间的计算，可存储起来供后续使用，
  # mergeCutHeight: 合并模块的阈值，越大模块越少,一般为0.25
  # minModuleSize: 每个模块里最少放多少个基因，设定越大模块越少
  # 输出结果根据模块中基因数目的多少，降序排列，依次编号为 `1-最大模块数`。
  # **0 (grey)**表示**未**分入任何模块的基因。
}

## 模块可视化，层级聚类树展示各个模块
# 模块构建完成后，将模块颜色labels转换为可视化所需的颜色，最后用层级聚类树展示各个模块
if(T){
  # Convert labels to colors for plotting
  moduleColors <- labels2colors(net$colors)
  table(moduleColors)
  # Plot the dendrogram and the module colors underneath
  pdf("step3_genes-modules_ClusterDendrogram.pdf",width = 16,height = 12)
  plotDendroAndColors(net$dendrograms[[1]], moduleColors[net$blockGenes[[1]]],
                      "Module colors",
                      dendroLabels = FALSE, hang = 0.03,
                      addGuide = TRUE, guideHang = 0.05)
  dev.off()
}
# 最后将得到的模块信息和模块颜色保存起来，方便后续的分析和可视化使用
save(net, moduleColors, file = "step3_genes_modules.Rdata")

这里面需要注意一个同名函数的报错

WGCNA包的blockwiseModules函数debug

https://www.jianshu.com/p/5b124641e4bc

我的分析结果中每个模块中基因的数量：

感觉作者module含的genes也不多：

5.关联基因模块与表型

终于要到那个有问题的相关性热图了

####################### 4.关联基因模块与表型 #####################################
rm(list = ls())  
load(file = "step1_input.Rdata")
load(file = "step2_power_value.Rdata")
load(file = "step3_genes_modules.Rdata")

##### 1.模块与表型的相关性热图----------------
if(T){
  datTraits$exp_group <- as.factor(datTraits$exp_group)
  # 表达矩阵对应的表型变量为设计矩阵，并设置组别
  design <- model.matrix(~0+datTraits$exp_group)
  colnames(design) <- levels(datTraits$exp_group) #获取组别信息
  # 计算所有模块的特征基因eigengenes.
  MES0 <- moduleEigengenes(datExpr,moduleColors)$eigengenes
  # MEblue      MEbrown       MEgrey   MEturquoise     MEyellow
  # LH_41_hg19 -0.054224196 -0.113190330  0.017756936 -0.0596382729  0.202767908
  # LH_42_hg19 -0.088915385 -0.133718403 -0.078586736 -0.0917636118 -0.026137421
  MEs <- orderMEs(MES0)  #将相关的特征基因放在一起
  # > dim(MEs)
  # [1] 40  5
  # > MEs    # 其实就换了下module的顺序
  # MEblue     MEyellow      MEbrown   MEturquoise       MEgrey
  # LH_41_hg19 -0.054224196  0.202767908 -0.113190330 -0.0596382729  0.017756936
  # LH_42_hg19 -0.088915385 -0.026137421 -0.133718403 -0.0917636118 -0.078586736
  
  # 可操作地方！！！！！
  # 计算模块与表型之间的相关性矩阵
  moduleTraitCor <- cor(MEs,design,use = "p")
  # > moduleTraitCor
  # [,1]       [,2]
  # MEblue       0.6998699 -0.6998699
  # MEyellow    -0.4405819  0.4405819
  # MEbrown      0.4076938 -0.4076938
  # MEturquoise  0.5932559 -0.5932559
  # MEgrey       0.1919625 -0.1919625
  moduleTraitPvalue <- corPvalueStudent(moduleTraitCor,nSamples)
  # [,1]         [,2]
  # MEblue      5.014994e-07 5.014994e-07
  # MEyellow    4.437355e-03 4.437355e-03
  # MEbrown     9.020902e-03 9.020902e-03
  # MEturquoise 5.474019e-05 5.474019e-05
  # MEgrey      2.353641e-01 2.353641e-01
  
  # 设定热图的文本标签为矩阵中每个元素的相关性值及其对应的显著性水平
  textMatrix <- paste0(signif(moduleTraitCor,2),"\n(",
                       signif(moduleTraitPvalue,1),")")
  # > textMatrix
  # [1] "0.7\n(5e-07)"   "-0.44\n(0.004)" "0.41\n(0.009)"  "0.59\n(5e-05)" 
  # [5] "0.19\n(0.2)"    "-0.7\n(5e-07)"  "0.44\n(0.004)"  "-0.41\n(0.009)"
  # [9] "-0.59\n(5e-05)" "-0.19\n(0.2)" 
  # signif 保留位数
  # dim函数用于更改或显示矩阵的维数，包括行数和列数。
  # 在这里，使用dim(textMatrix)函数来设置textMatrix的维度，
  # 其行数和列数与模块与表型相关性的矩阵相同，以确保文本正确地添加到每个网格中。
  dim(textMatrix) <- dim(moduleTraitCor)
  # [,1]             [,2]            
  # [1,] "0.7\n(5e-07)"   "-0.7\n(5e-07)" 
  # [2,] "-0.44\n(0.004)" "0.44\n(0.004)" 
  # [3,] "0.41\n(0.009)"  "-0.41\n(0.009)"
  # [4,] "0.59\n(5e-05)"  "-0.59\n(5e-05)"
  # [5,] "0.19\n(0.2)"    "-0.19\n(0.2)" 
  
  pdf("step4_Module-trait-relationship_heatmap.pdf",
      width = 2*length(colnames(design)), 
      height = 0.6*length(names(MEs)) )
  par(mar=c(5, 9, 3, 3)) #留白：下、左、上、右
  # 绘制
  labeledHeatmap(Matrix = moduleTraitCor,  # 要绘制热图的矩阵
                 xLabels = colnames(design),  # X轴上的标签
                 yLabels = names(MEs),  # Y轴上的标签
                 ySymbols = names(MEs),  # 线条的标记符号
                 colorLabels = F,  # 颜色标签 如果为TRUE，则在热图下方添加一个注释区，表示颜色对应的值的范围
                 colors = blueWhiteRed(50),  # 颜色网格的颜色图谱
                 textMatrix = textMatrix,  # 添加到每个网格中的文本（可以是数字或字符）
                 setStdMargins = F,  # 是否将margins调整为标准大小
                 cex.text = 0.5,  # 文本大小
                 zlim = c(-1,1),   # 颜色图谱的值范围
                 main = "Module-trait relationships")  # 热图的标题
  dev.off()
  save(design, file = "step4_design.Rdata")
}

这个是没有加入其他表型信息的相关性热图：

加入样本其它表型信息

样本信息的加入也是通过相关性系数矩阵来的

手动设计前面的二分组实验矩阵进行相关性计算也一样

我们就将以此为基础，往里加其它表型-模块的相关性系数

获取其它表型信息：

试下getGEO能不能直接获取到，不行：

后在数据集文章网页找到：

https://www.ahajournals.org/doi/10.1161/CIRCRESAHA.117.311572

缺血性中风是指脑血管发生堵塞或狭窄，导致脑部缺血（血液供应不足）的情况。“Infraction volume”（梗死体积）是指中风后形成的脑部组织梗死区域的大小。

我本来以为这些表型绝大多数都是数值型变量，应该是拿连续性变量和module进行相关性分析，但我细细一看，好家伙，这些表型信息作者用的都是根据原有二分组得到的平均或者百分比结果！那这些表型根本没反映对应的样品各自具体的表型信息啊，全都反映的是原有的实验分组后的平均表型。。。。

可这样所有这些表型还有意义吗，不都是根据IS vs HC 分组的吗

哦！我可算知道你这个图为啥这么奇怪了

之所以，所有的表型在module上的相关性趋势都是一样的，并且大部分相关性系数大小一样（除了sex和后面计算的这个分数，具体如何改变的因为作者没有透露，我们无从知晓，但趋势仍一样），是因为实际上都是根据IS vs HC的分组！

黄色框起来的部分和没框起来的刚好相关性正负相反，如果你把这些表型，如sex对应的male、female中选择的背景颠倒（稀疏矩阵中 0 1交换），就会得到module趋势完全一样的相关性

可以看到除了性别这个分类变量（作者应该是手动把Female改成了sex），其余黄色框起来的和没框起来的数值型变量大小关系和作者这张模块与表型相关性热图一致

就像前面，可以看到HC IS对应的相关性系数是大小相同正负相反，因为他们是一组对照

我的结果：

# 数据集文章有

if(T){
  # 计算所有模块的特征基因eigengenes.
  MES0 <- moduleEigengenes(datExpr,moduleColors)$eigengenes
  # MEblue      MEbrown       MEgrey   MEturquoise     MEyellow
  # LH_41_hg19 -0.054224196 -0.113190330  0.017756936 -0.0596382729  0.202767908
  # LH_42_hg19 -0.088915385 -0.133718403 -0.078586736 -0.0917636118 -0.026137421
  MEs <- orderMEs(MES0)  #将相关的特征基因放在一起
  # > dim(MEs)
  # [1] 40  5
  # > MEs    # 其实就换了下module的顺序
  # MEblue     MEyellow      MEbrown   MEturquoise       MEgrey
  # LH_41_hg19 -0.054224196  0.202767908 -0.113190330 -0.0596382729  0.017756936
  # LH_42_hg19 -0.088915385 -0.026137421 -0.133718403 -0.0917636118 -0.078586736
  
  # 可操作地方！！！！！
  # 计算模块与表型之间的相关性矩阵
  phenotyes <- read.table("phenotypes.csv",header = TRUE,sep = ",")
  selected <- c(1,2,3,4,5,6,7,9,16,19,32)
  phenotyes <- phenotyes[selected,]
  phenotyes$HC[11] <- 0
  rownames(phenotyes) <- phenotyes$X
  phenotyes$X <- NULL
  phenotyes[2,] <- rev(phenotyes[2,])
  rn <- rownames(phenotyes)
  rn[2] <- "sex"
  rownames(phenotyes) <- rn
  
  all_cor_textMat <- apply(phenotyes,1,function(x){
    design <- rep(x,each=20)
    colnames(design) <- rownames(x) #获取组别信息
    # 要求 equal numbers of rows 行数一致 design 40 * 1
    moduleTraitCor <- cor(MEs,design,use = "p")  # 5module * 1phenotype
    moduleTraitPvalue <- corPvalueStudent(moduleTraitCor,nSamples)
    textMatrix <- paste0(signif(moduleTraitCor,2),"\n(",
                         signif(moduleTraitPvalue,1),")")
    # dim(textMatrix) <- dim(moduleTraitCor)
    return(textMatrix)
  })
  
  
  all_cor_moduleTrtCor <- apply(phenotyes,1,function(x){
    design <- rep(x,each=20)
    colnames(design) <- rownames(x) #获取组别信息
    # 要求 equal numbers of rows 行数一致 design 40 * 1
    moduleTraitCor <- cor(MEs,design,use = "p")  # 5module * 1phenotype
    return(moduleTraitCor)
  })

  
  pdf("step4_ALLPhenotype_Module-relationship_heatmap.pdf",
      width = 1.4*length(rownames(phenotyes)), 
      height = 0.6*length(names(MEs)) )
  par(mar=c(5, 9, 3, 3)) #留白：下、左、上、右
  # 绘制
  labeledHeatmap(Matrix = all_cor_moduleTrtCor,  # 要绘制热图的矩阵
                 xLabels = rownames(phenotyes),  # X轴上的标签
                 yLabels = names(MEs),  # Y轴上的标签
                 ySymbols = names(MEs),  # 线条的标记符号
                 colorLabels = F,  # 颜色标签 如果为TRUE，则在热图下方添加一个注释区，表示颜色对应的值的范围
                 colors = blueWhiteRed(50),  # 颜色网格的颜色图谱
                 textMatrix = all_cor_textMat,  # 添加到每个网格中的文本（可以是数字或字符）
                 setStdMargins = F,  # 是否将margins调整为标准大小
                 cex.text = 0.5,  # 文本大小
                 zlim = c(-1,1),   # 颜色图谱的值范围
                 main = "Module-trait relationships")  # 热图的标题
  dev.off()
}