R语言分析3：加权相关网络分析(WGCNA)

lichengxin • • 杂文

加权基因共表达网络分析（WGCNA, Weighted gene co-expression network analysis）是一种用来描述不同样品之间基因关联模式的系统生物学方法，可以用来鉴定高度协同变化的基因集, 并根据基因集的内连性和基因集与表型之间的关联鉴定候补生物标记基因或治疗靶点。相比于只关注差异表达的基因，WGCNA利用数千或近万个变化最大的基因或全部基因的信息识别感兴趣的基因集，并与表型进行显著性关联分析。一是充分利用了信息，二是把数千个基因与表型的关联转换为数个基因集与表型的关联，免去了多重假设检验校正的问题。

专有名词：
ꔷ 共表达网络：定义为加权基因网络。点代表基因，边代表基因表达相关性
ꔷ Module(模块)：高度內连的基因集。在无向网络中，模块内是高度相关的基因；在有向网络中，模块内是高度正相关的基因。把基因聚类成模块后，可以对每个模块进行三个层次的分析：1. 功能富集分析查看其功能特征是否与研究目的相符；2. 模块与性状进行关联分析，找出与关注性状相关度最高的模块；3. 模块与样本进行关联分析，找到样品特异高表达的模块
ꔷ 连接度(Connectivity)：类似于网络中 "度"(degree)的概念。每个基因的连接度是与其相连的基因的边属性之和
ꔷ 邻接矩阵(Adjacency matrix)：基因和基因之间的加权相关性值构成的矩阵
ꔷ Intramodular connectivity：模内连通性测量，给定基因相对于特定模块的基因如何连接或共表达，判断基因所属关系
ꔷ Module eigengene E：给定模型的第一主成分，代表整个模型的基因表达谱，可很好的用一个向量代替了一个矩阵，方便后期计算
ꔷ Eigengene significance：当获得微阵列样品特征y（例如病例对照状态或体重）时，可以将模块特征基因与此结果相关联， 相关系数称为特征基因显著性
ꔷ Module Membership KME ：对于每个基因，通过将其基因表达谱与给定模块的模块本征基因相关联来定义模块成员的“模糊”度量。如，MM^{blue}(i) = K_{cor,r}^{blue} = cor(x^i,E^{blue})用来测量基因i与蓝色模块特征基因的相关程度。若MM^{blue}(i)接近0，则第i个基因不属于blue模块的一部分，相反，MM^{blue}(i)接近于1或 -1，则它与蓝色模块基因高度相关（正或负）。
ꔷ Gene significance GS：GS_i的绝对值越高，第i个基因的生物学意义就越大。
ꔷ Module significance：给定模块中所有基因的平均绝对基因显著性的度量。 当将基因显着性定义为基因表达与外部性状y的相关性时，此度量往往与模块特征基因与y的相关性高度相关
ꔷ TOM (Topological overlap matrix)：把邻接矩阵转换为拓扑重叠矩阵，以降低噪音和假相关，获得的新距离矩阵，这个信息可拿来构建网络或绘制TOM图
ꔷ Hub gene：关键基因 (连接度最多或连接多个模块的基因)
ꔷ 软阈值：相关性值进行幂次运算幂次的值也就是软阈值;

WGCNA-1

1. 数据输入、清洗和预处理

1.1 加载基因表达数据

library(WGCNA)

dim(FPKM) # 仍选取TCGA- LUAD的FPKM数据
## [1] 59427   600

## 筛选中位绝对偏差前75%的基因，至少MAD大于0.01
## 筛选后会降低运算量，也会失去部分信息
## 也可不做筛选，使MAD大于0即可
m.mad <- apply(FPKM, 1, mad)
dataExprVar <- FPKM[which(m.mad > max(quantile(m.mad, probs = seq(0, 1, 0.25))[2], 0.01)),]
# 过滤基因，选取绝对中位差Top5000的基因
# dataExprVar <- t(exp[order(apply(FPKM, 1, mad), decreasing = T)[1:5000],])
dim(dataExprVar)
## [1] 30976   600

## 转换为样品在行，基因在列的矩阵
LUAD_Expr0 <- as.data.frame(t(dataExprVar))

★ 注意：并不推荐使用差异基因作为输入矩阵，通过差异表达基因过滤将会导致一个（或几个高度相关的）基因聚成一个模块，同时，也破坏了无标度拓扑的假设，所以通过无标度拓扑拟合来选择软阈值的将会失败

1.2 数据清洗

## 检测缺失值
gsg = goodSamplesGenes(LUAD_Expr0, verbose = 3)
gsg$allOK;
## 1] TRUE

# 若样本检查语句返回TRUE，那么没有缺失值。否则，可通过以下代码来移除缺失值
if (!gsg$allOK){
  # Optionally, print the gene and sample names that were removed:
  if (sum(!gsg$goodGenes)>0) 
    printFlush(paste("Removing genes:", 
                     paste(names(LUAD_Expr0)[!gsg$goodGenes], collapse = ",")));
  if (sum(!gsg$goodSamples)>0) 
    printFlush(paste("Removing samples:", 
                     paste(rownames(LUAD_Expr0)[!gsg$goodSamples], collapse = ",")));
  # Remove the offending genes and samples from the data:
  LUAD_Expr0 = LUAD_Expr0[gsg$goodSamples, gsg$goodGenes]
}

## 对样本进行聚类（与随后的基因聚类相比），看看是否有明显的异常值
# hclusts聚类算法, dist计算基因之间的距离
sampleTree = hclust(dist(LUAD_Expr0), method = "average");
# 设置文字大小
par(cex = 0.2);
# 设置图像边距c(bottom, left, top, right) 
par(mar = c(0,4,2,0))
# 画图 main标题，sub子标题，xlab x轴标题，cex.lab标题字体大小，cex.axis坐标轴刻度大小，cex.main主标题字体
plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5, cex.axis = 1.5, cex.main = 5)

# 在上图上画红线
abline(h = 1e+05, col = "red");
# Determine cluster under the line
# 剪枝算法，cutHeight 修剪树枝的高度 minSize集群最小数
clust = cutreeStatic(sampleTree, cutHeight = 1e+05, minSize = 5)
# 剪枝结果
table(clust)
## clust
##   0   1   2 
##   2 591   7 

# clust 1 contains the samples we want to keep
keepSamples = (clust==1)
# 符合要求的数据
LUAD_Expr = LUAD_Expr0[keepSamples, ]

# 提取行和列
nSamples = nrow(LUAD_Expr)
nGenes = ncol(LUAD_Expr)
WGCNA-2

★ 此处去除左侧9个样本

1.3 加载临床特征数据

# 读取临床数据
LUAD_clin0 <- jsonlite::fromJSON("/data/shumin/Cibersort/LUAD/clinical.cart.2023-06-08.json")
n = length(LUAD_clin0$diagnoses)

# 提取数据
case_id = classfication_of_tumor = c(rep(0, n))
tumor_stage = gender = c(rep(0, n))
year_to_birth = year_to_death =  c(rep(0, n))
year_to_diagnosis = days_to_death = c(rep(0, n))
age = deadORlive = race = alcohol = smoked = c(rep(0, n))

for (i in 1:n) {
  case_id[i] = LUAD_clin0$case_id[[I]]
  classfication_of_tumor[i] = LUAD_clin0$diagnoses[[i]]$classification_of_tumor
  tumor_stage[i] = LUAD_clin0$diagnoses[[i]]$tumor_grade
  gender[i] = LUAD_clin0$demographic$gender[[I]]
  year_to_birth[i] = ifelse(
    is.null(LUAD_clin0$demographic$year_of_birth[[i]]),
    "notReport",
    LUAD_clin0$demographic$year_of_birth[[I]]
  )
  year_to_death[i] = ifelse(
    is.null(LUAD_clin0$demographic$year_of_death[[i]]),
    "notReport",
    LUAD_clin0$demographic$year_of_death[[I]]
  )
  year_to_diagnosis[i] = ifelse(
    is.null(LUAD_clin0$diagnoses[[i]]$year_of_diagnosis),
    "notReport",
    LUAD_clin0$diagnoses[[i]]$year_of_diagnosis
  )
  days_to_death[i] = ifelse(
    is.null(LUAD_clin0$demographic$days_to_death[[i]]),
    "notReport",
    LUAD_clin0$demographic$days_to_death[[I]]
  )
  age[i] = ifelse(
    is.null(LUAD_clin0$demographic$age_at_index[[i]]),
    "notReport",
    LUAD_clin0$demographic$age_at_index[[I]]
  )
  deadORlive[i] = ifelse(
    is.null(LUAD_clin0$demographic$vital_status[[i]]),
    "notReport",
    LUAD_clin0$demographic$vital_status[[I]]
  )
  race[i] = ifelse(
    is.null(LUAD_clin0$demographic$race[[i]]),
    "notReport",
    LUAD_clin0$demographic$race[[I]]
  )
  alcohol[i] = ifelse(
    is.null(LUAD_clin0$exposures[[i]]$alcohol_history),
    "notReprot",
    LUAD_clin0$exposures[[i]]$alcohol_history
  )
  smoked[i] = ifelse(
    is.null(LUAD_clin0$exposures[[i]]$years_smoked),
    "notReport",
    LUAD_clin0$exposures[[i]]$years_smoked
  )
}

LUAD_clin1 <- data.frame(
  case_id,
  classfication_of_tumor,
  tumor_stage,
  gender,
  year_to_birth,
  year_to_death,
  year_to_diagnosis,
  days_to_death,
  age,
  deadORlive,
  race,
  alcohol,
  smoked
)

head(LUAD_clin1)
##                                case_id classfication_of_tumor  tumor_stage gender year_to_birth year_to_death year_to_diagnosis days_to_death age deadORlive                      race      alcohol smoked
##  1 a3fd20b2-e001-44ab-9716-754e5ae70808           not reported not reported female          1935            NA              2012            NA  77      Alive                     white Not Reported     NA
##  2 25d4ea9e-f773-4f11-bac9-64efdad73211           not reported not reported female          1950            NA              2009            NA  59      Alive                     white Not Reported     45
##  3 a52e99d6-a61a-439d-b0b1-ca7a0eabcb04           not reported not reported female          1969            NA              2011            NA  42      Alive                     white Not Reported      2
##  4 27fceec1-3298-4cdd-a4e6-8f5cf34604f0           not reported not reported   male          1957            NA              2010            NA  53      Alive black or african american Not Reported     NA
##  5 2923e404-38f2-437a-b57e-23401fbe0273           not reported not reported   male          1965            NA              2011            NA  46      Alive                     white Not Reported     NA
##  6 a65700c2-e58c-4fd4-aeb1-5686b8f4d212           not reported not reported   male          1935            NA              2009            NA  74      Alive                     white Not Reported     NA

# 读取meta数据（用于id转换）
LUAD_meta <- jsonlite::fromJSON("/data/shumin/Cibersort/LUAD/metadata.cart.2023-06-08.json")
sample_id <- sapply(LUAD_meta$associated_entities, function(x){x[,1]})
case_id <- sapply(LUAD_meta$associated_entities, function(x){x[,3]})
anno_data <- data.frame(sample_id, case_id) 

head(anno_data)
##                      sample_id                              case_id
##  1 TCGA-44-8120-01A-11R-2241-07 83e38dbd-edab-47f2-b19f-6ea38fc6bece
##  2 TCGA-99-8025-01A-11R-2241-07 84c3ba70-afa7-4b69-be69-7ec8d6022c56
##  3 TCGA-50-6594-01A-11R-1755-07 8504fd86-a70a-4cba-9ec8-25c9e60ca549
##  4 TCGA-78-7161-01A-11R-2039-07 86faf16c-56fd-4b7b-a6b2-b4c83bec93e1
##  5 TCGA-55-7903-01A-11R-2170-07 77c4dbb2-eceb-4e0d-bcde-63dc817d5f35
##  6 TCGA-38-4632-01A-01R-1755-07 875333ab-9048-462d-aaa2-693ad127e3cc

LUAD_clin2 <- merge(anno_data, LUAD_clin1, by = "case_id")

LUAD_clin3 <- LUAD_clin2[, c("sample_id","gender","year_to_birth","year_to_death","year_to_diagnosis","days_to_death", "age","deadORlive","race","smoked")] # 选取后续分析需要的列
LUAD_clin3$gender <-  ifelse(LUAD_clin3$gender == 'female', 1, 2)
LUAD_clin3$deadORlive <- ifelse(LUAD_clin3$deadORlive == 'Alive', 0, 1)
LUAD_clin3$race <- ifelse(LUAD_clin3$race == 'not reported', 0, 
                          ifelse (LUAD_clin3$race == 'american indian or alaska native', 1,
                                  ifelse(LUAD_clin3$race == 'asian', 2,
                                         ifelse(LUAD_clin3$race == 'black or african american', 3, 4))))

# 形成一个类似于表达数据的数据框架，以保存临床特征
# 提取行名
LUAD_Samples = rownames(LUAD_Expr)
# 数据匹配 返回匹配行
Trait_Rows = match(LUAD_Samples, LUAD_clin3$sample_id);
# 提取指定要求行
data_Traits = LUAD_clin3[Trait_Rows, -1];
# 提取行名
rownames(data_Traits) = LUAD_clin3[Trait_Rows, 1];
# 垃圾回收
collectGarbage();

# Re-cluster samples
# 画聚类图
LUAD_Tree = hclust(dist(LUAD_Expr), method = "average")
# Convert traits to a color representation: white means low, red means high, grey means missing entry
# 画表型的热图
# 将特征转换为颜色表示：白色表示低，红色表示高，灰色表示缺少条目
# 如果signed为true 以绿色开头代表最大负值，以白色开头代表零附近的值，然后变为红色代表正值
traitColors = numbers2colors(data_Traits, signed = FALSE);
# Plot the sample dendrogram and the colors underneath.
# 绘制出树状图和下面的颜色 
plotDendroAndColors(LUAD_Tree, traitColors, groupLabels = names(data_Traits), dendroLabels = FALSE, main = "LUAD dendrogram and trait heatmap")
WGCNA-3

2. 建设表达网络与模块检测

此步骤是使用WGCNA方法进行所有网络分析的基础。WGCNA提出三种不同的方法构建网络并识别模块：

/ ꔷ 使用方便的一步网络结构和模块检测功能，适合希望以最小努力达到结果的用户；
/ ꔷ 为希望使用定制/替代方法进行实验的用户逐步构建网络和模块检测；
/ ꔷ 一种自动分块网络结构和模块检测方法，适用于希望分析太大而无法同时分析的数据集的用户。

2.1 自动一步构建网络与模块检测

2.1.1 软阈值的选择：网络拓扑分析

构建一个加权基因网络需要选择软阈值幂β来计算邻接矩阵权重参数，即将基因间的相关系数进行乘方运算来表征其相关性，首先需要确定乘方的值

# Choose a set of soft-thresholding powers
# 给出候选的β值，c(1:10)表示1到10；seq(from = 12, to = 20, by = 2)表示从12开始间隔两个数到20
powers = c(c(1:10), seq(from = 12, to = 20, by = 2))
powers
# Call the network topology analysis function 调用网络拓扑分析函数
# verbose表示输出结果详细程度
sft = pickSoftThreshold(LUAD_Expr, powerVector = powers, verbose = 0);

# sft这中保存了每个powers值计算出来的网络特征,其中powerEstimate就是最佳power值，fitIndices保存了每个power对应的网络的特征。
str(sft)
## List of 2
##  $ powerEstimate: num 2
##  $ fitIndices   :'data.frame':   15 obs. of  7 variables:
##   ..$ Power         : num [1:15] 1 2 3 4 5 6 7 8 9 10 ...
##   ..$ SFT.R.sq      : num [1:15] 0.123 0.916 0.916 0.932 0.938 ...
##   ..$ slope         : num [1:15] -0.295 -1 -0.996 -0.964 -0.954 ...
##   ..$ truncated.R.sq: num [1:15] 0.393 0.921 0.94 0.933 0.925 ...
##   ..$ mean.k.       : num [1:15] 4171 1288 613 364 244 ...
##   ..$ median.k.     : num [1:15] 3579.7 676.1 173.7 54 19.4 ...
##   ..$ max.k.        : num [1:15] 9260 5325 3783 2938 2403 ...

# Plot the results 结果绘图
# 设置图的显示一行两列
par(mfrow = c(1,2));
cex1 = 0.9;
# 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.90, col = "red")

# sft$fitIndices 保存了每个power构建的相关性网络中的连接度的统计值，k就是连接度值，每个power值提供了max, median, max3种连接度的统计量
# 对连接度的均值进行可视化
# 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")

WGCNA-4

★ R^{2}为无尺度网络评价系数，一般设置为0.9， β取值标准：R^{2}第一次到达0.9时对应的β值
★ 横坐标为软阈值的梯度，第一幅图的纵坐标为无标度网络适应系数，越大越好；第二幅图的纵坐标为节点的平均连通度，越小越好

2.1.2 一步构建网络与模块检测

# 打开多线程
enableWGCNAThreads(nThreads = 10) 

# LUAD_Expr表达数据，TOMType拓扑重叠矩阵计算方式，minModuleSize用于模块检测的最小模块尺寸,
# reassignThreshold 是否在模块之间重新分配基因的p值比率阈值，mergeCutHeight 树状图切割高度
# numericLabels 返回的模块应该用颜色（FALSE）还是数字（TRUE）标记,pamRespectsDendro树状图相关参数
# saveTOMs 字符串的向量，saveTOMFileBase 包含包含共识拓扑重叠文件的文件名库的字符串
net = blockwiseModules(LUAD_Expr, power = sft$powerEstimate, 
                       maxBlockSize = 2000, # 最大模块数量
                       TOMType = "unsigned", minModuleSize = 30, reassignThreshold = 0, # 构建无尺度网络，最小模块基因数为30
                       mergeCutHeight = 0.25, # 需要合并模块的阈值
                       numericLabels = TRUE, # 以数字作为模块的名字
                       pamRespectsDendro = FALSE, saveTOMs = TRUE,
                       saveTOMFileBase = "LUADTOM", verbose = 3)

# 回到网络分析，查看标识了多少个模块以及模块大小
table(net$colors)
## 
##    0    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16   17   18   19   20   21   22   23   24   25   26   27   28   29   30   31   32   33   34   35   36   37   38   39   40   41   42   43   44   45   46   47   48   49   50   51   52   53   54   55   56 
##  4348 9126 1444 1202 1158 1141  969  725  581  531  444  425  391  359  324  317  280  274  263  252  225  186  178  152  150  142  137  128  125  118  115  115  114  112  112  108  108  106  100   96   96   92   92   90   89   88   83   82   81   81   81   80   80   79   78   78   77 
##    57   58   59   60   61   62   63   64   65   66   67   68   69   70   71   72   73   74   75   76   77   78   79   80   81   82   83   84   85   86   87   88   89   90   91   92   93   94   95   96   97   98   99  100  101  102  103  104 
##    74   74   74   68   67   64   64   64   62   61   60   59   58   57   56   55   55   54   52   52   52   52   50   50   50   49   49   49   47   45   43   43   39   39   38   38   38   36   36   36   34   34   33   33   33   31   31   30 

# 将标签转化为绘图颜色
moduleColors = labels2colors(net$colors)
# 绘制树状图和下面的模块颜色
# dendroLabels树状图标签。设置为FALSE完全禁用树状图标签；设置为NULL使用的行标签LUAD_Expr
# addGuide是否应在树状图中添加垂直的“指导线”？线条使识别单个样本的颜色代码更加容易。
plotDendroAndColors(net$dendrograms[[1]], moduleColors[net$blockGenes[[1]]], "Module colors",
                    dendroLabels = FALSE, hang = 0.03, addGuide = TRUE, guideHang = 0.05)
WGCNA-5

★ 指示有104个模块，按大小降序标记为1至104，大小范围为9126至30个基因； 标签0保留用于所有模块外部的基因
★ 无法聚类到模块中的基因会标示为灰色，如果灰色区域较多，可能由于样本中基因共表达趋势不明显，可能需要调整基因过滤的方法

# 保存后续分析所需的模块分配和模块本征信息，可由recutBlockwiseTrees函数应用修改后的条件而不必重新计算网络和聚类树状图
moduleLabels = net$colors
moduleColors = labels2colors(net$colors)
MEs = net$MEs;
geneTree = net$dendrograms[[1]];
save(MEs, moduleLabels, moduleColors, geneTree, file = "TCGA-LUAD-02-networkConstruction-auto.RData")

2.2 分步网络构建和模块检测

https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/Tutorials/FemaleLiver-02-networkConstr-man.pdf

2.3 逐块网络构建和模块检测

https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/Tutorials/FemaleLiver-02-networkConstr-blockwise.pdf

3. 分步网络构建和模块检测

3.1 量化模块-特质关联

在此分析中，我们想确定与所测量的临床特征显着相关的模块。 由于我们已经为每个模块建立了一个概要文件（特征基因），因此我们只需将特征基因与外部特征相关联，然后寻找最重要的关联，实际上是计算模块的ME值与表型的相关系数：

# Define numbers of genes and samples
# 获得基因数和样本数
nGenes = ncol(LUAD_Expr);
nSamples = nrow(LUAD_Expr);

# Recalculate MEs with color labels
# 用彩色标签重新计算MEs
# 在给定的单个数据集中计算模块的模块本征基因
MEs0 = moduleEigengenes(LUAD_Expr, moduleColors)$eigengenes
# 对给定的（特征）向量进行重新排序，以使相似的向量（通过相关性度量）彼此相邻
MEs = orderMEs(MEs0)

# 计算module的ME值与表型的相关系数（pearson）
moduleTraitCor = cor(MEs, data_Traits, use = "p");
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples);
# 使用的是其他相关性方法，则可以使用 bicorAndPvalue 函数来计算显著性
# modTraitCorP = bicorAndPvalue(MEs_col, design)
# modTraitCor = modTraitCorP$bicor
# modTraitP   = modTraitCorP$p

names(MEs)
##   [1] "MEdarkgrey"        "MEyellowgreen"     "MEivory"           "MEnavajowhite1"    "MEindianred3"      "MElightcyan"       "MEgreenyellow"     "MEskyblue1"        "MEbisque4"         "MEblack"           "MEbrown"           "MEmagenta"         "MEmidnightblue"    "MEpalevioletred3" 
##   [15] "MEskyblue4"        "MEmediumpurple1"   "MEthistle1"        "MEblue2"           "MEorangered3"      "MEdarkviolet"      "MEplum1"           "MEthistle"         "MEdarkred"         "MEwhite"           "MEfirebrick3"      "MEblueviolet"      "MEsienna4"         "MElightsteelblue" 
##   [29] "MElightyellow"     "MEpalevioletred2"  "MEcoral3"          "MEpink4"           "MEhoneydew"        "MEyellow3"         "MEnavajowhite2"    "MEcoral1"          "MEtan"             "MEblue"            "MEdarkolivegreen"  "MEthistle3"        "MEdarkseagreen4"   "MEplum3"          
##   [43] "MEbrown2"          "MElavenderblush3"  "MEhoneydew1"       "MElightslateblue"  "MEdarkolivegreen4" "MEdarkturquoise"   "MEviolet"          "MEcoral"           "MEpink"            "MEmagenta4"        "MEfirebrick4"      "MEpurple"          "MEred"             "MElightcoral"     
##   [57] "MEyellow"          "MEmediumpurple4"   "MEskyblue"         "MEcyan"            "MEdarkgreen"       "MElightpink4"      "MEskyblue2"        "MEmediumorchid"    "MElightsteelblue1" "MEmaroon"          "MEmediumpurple2"   "MEplum2"           "MEsteelblue"       "MEdarkorange2"    
##   [71] "MEplum"            "MEdarkmagenta"     "MEorange"          "MEdarkorange"      "MEsalmon"          "MElightgreen"      "MEturquoise"       "MEroyalblue"       "MEyellow4"         "MEorangered4"      "MEpaleturquoise"   "MElightblue4"      "MEindianred4"      "MEbrown4"         
##   [85] "MEsalmon4"         "MEsienna3"         "MEthistle2"        "MEgrey60"          "MEgreen"           "MEorangered1"      "MEmediumpurple3"   "MEdarkslateblue"   "MElightpink3"      "MEantiquewhite2"   "MElightcyan1"      "MEfloralwhite"     "MElavenderblush2"  "MEdarkolivegreen2"
##   [99] "MEskyblue3"        "MEsalmon2"         "MEcoral2"          "MEsaddlebrown"     "MEantiquewhite4"   "MEdarkseagreen3"   "MEgrey" 

# sizeGrWindow(10,6)
# 显示相关性及其p值
textMatrix = paste(signif(moduleTraitCor, 2), "/n(",signif(moduleTraitPvalue, 1), ")", sep = "");
dim(textMatrix) = dim(moduleTraitCor)
# Display the correlation values within a heatmap plot/
# ySymbols 当ylabels使用时所使用的其他标签； colorLabels 应该使用颜色标签吗
# colors 颜色； textMatrix 单元格名字
labeledHeatmap(Matrix = moduleTraitCor, xLabels = names(data_Traits), yLabels = names(MEs), ySymbols = names(MEs),
               colorLabels = FALSE,colors = greenWhiteRed(50), textMatrix = textMatrix,setStdMargins = FALSE,
               cex.text = 0.2, zlim = c(-1,1),
               main = paste("Module-trait relationships"))

WGCNA-6

★ 计算出模块与表型之间相关性之后，可以挑选最相关的那些模块来进行后续分析。但是，模块本身可能还包含很多的基因，还需要进一步识别关键基因

3.2 基因与性状和重要模块的关系：基因重要性和模块成员

量化阵列上所有基因与每个模块的相似性寻找重要模块:

# 定义包含数据特征权重列的变量权重
days_to_death = as.data.frame(data_Traits$days_to_death);
names(days_to_death) = "days_to_death"
geneModuleMembership = as.data.frame(cor(LUAD_Expr, MEs, use = "p"));
# 模块的名称(颜色)substring提取文本从第3个字母开始（此处输入法有问题，不用"#")
# modNames = substring(names(MEs), 3)
# 基因和模块的相关系数
geneModuleMembership = as.data.frame(cor(LUAD_Expr, 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和性状的关系
geneTraitSignificance = as.data.frame(cor(LUAD_Expr, days_to_death, use = "p"));
GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples));
names(geneTraitSignificance) = paste("GS.", names(days_to_death), sep="");
names(GSPvalue) = paste("p.GS.", names(days_to_death), sep = "");

3.3 模内分析：鉴定具有高GS和MM的基因

使用GS和MM度量，我们可以鉴定出对“days_to_death”以及在感兴趣的模块中具有较高模块成员性具有重要意义的基因。 例如，我们看一下与“days_to_death”关联“brown”模块。 我们在“brown”模块中绘制了基因重要性与模块成员关系的散点图。 在此模块中，GS和MM之间存在高度显着的负相关性

# 模型颜色
module = "brown"
# 匹配列
column = match(module, modNames);
moduleGenes = moduleColors==module;
#sizeGrWindow(7, 7);
par(mfrow = c(1,1));
# MM <- abs(geneModuleMembership[moduleGenes, column])
# GS <- abs(geneTraitSignificance[moduleGenes, 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)

WGCNA-7

★ GS是Gene Signicance，描述的是某一个基因与性状的相关性
★ MM是Module Membership，描述的是某一个基因与模块之间的相关性

3.4 导出模块网络

file <- "~/TCGA-LUAD.net"
TOM <- TOMsimilarityFromExpr(LUAD_Expr, power = 2, networkType = "unsigned")
dimnames(TOM) <- list(colnames(LUAD_Expr), colnames(LUAD_Expr))
modTOM <- TOM[moduleGenes, moduleGenes]

cyt <- exportNetworkToCytoscape(
  modTOM,
  edgeFile = paste0(file, module, ".edges.txt"),
  nodeFile = paste0(file, module, ".nodes.txt"),
  weighted = TRUE,
  threshold = 0, # threshold 默认为0.5, 可以根据自己的需要调整
  nodeNames = moduleGenes,
  nodeAttr = module
)

4. 使用WGCNA进行网络可视化

4.1 显示基因网络

可视化加权网络的一种方法是绘制其热图，热图的每一行和每一列都对应一个基因。 热图可以描述邻接或拓扑重叠，浅色表示低邻接（重叠），而深色表示更高的邻接（重叠）。 另外，沿着热图的顶部和左侧绘制了基因树状图和模块颜色。

# 重新计算拓扑重叠：通过保存TOM可以更有效地完成此操作
dissTOM = 1-TOMsimilarityFromExpr(LUAD_Expr, power = 2);
# 变换dissTOM（对dissTOM矩阵进行指数转换，使中等强度的关系更容易在热图上展示出来）；
plotTOM = dissTOM^7;
# #将对角线数据设为NA
diag(plotTOM) = NA;
# Call the plot function
# sizeGrWindow(9,9)
# 基因的聚类树聚类时的距离为1-TOM值结合基因间的距离，即1-TOM值，用热图展示
# TOMplot(plotTOM, geneTree, moduleColors, main = "Network heatmap plot, all genes") # 此处我的基因数量过多？会报错：Error in x[, iy] : subscript out of bounds

★ 数据不合理，可根据探针集或基因的平均表达量或方差（如中位数或绝对中位差）重新进行过滤，低表达或无变化的基因通常代表噪音。

# 限制基因数量以加快绘图速度。 注意，一个基因子集的基因树状图通常看起来与所有基因的基因树状图不同：
nSelect = 400
# 为了反复重现结果, 这里设置了随机种子；
set.seed(10);
select = sample(nGenes, size = nSelect);
# 提取抽取基因相应的TOM矩阵
selectTOM = dissTOM[select, select];
# 重新画聚类图
selectTree = hclust(as.dist(selectTOM), method = "average")
selectColors = moduleColors[select];
# Open a graphical window
# sizeGrWindow(9,9)
# Taking the dissimilarity to a power, say 10, makes the plot more informative by effectively changing
# the color palette; setting the diagonal to NA also improves the clarity of the plot
plotDiss = selectTOM^7;
diag(plotDiss) = NA;
# 更换配色
# mycolor <- gplots::colorpanel(250, 'red', "orange", 'lemonchiffon')
TOMplot(plotDiss, selectTree, selectColors,
        main = "Network heatmap plot, selected genes",
        # col = mycolor)
WGCNA-8

4.2 可视化特征基因网络

研究找到的模块之间的关系通常很有趣。 可以使用特征基因作为代表特征，并通过特征基因相关性来量化模块相似性。 该软件包包含一个方便的函数plotEigengeneNetworks，该函数生成特征基因网络的摘要图。 通常，向特征基因添加临床特征（或多个特征）以了解特征如何适合特征基因网络是有益的

# Recalculate module eigengenes
# 重新计算基因特征值
MEs = moduleEigengenes(LUAD_Expr, moduleColors)$eigengenes
# Isolate weight from the clinical traits
days_to_death = as.data.frame(data_Traits$days_to_death);
names(days_to_death) = "days_to_death"
# Add the weight to existing module eigengenes
MET = orderMEs(cbind(MEs, days_to_death))
# Plot the relationships among the eigengenes and the trait
# sizeGrWindow(5,7.5);
par(cex = 0.2)
# 画树形图
# marDendro给出树状图的边距设置，marHeatmap热图边距设置
plotEigengeneNetworks(MET, "",  marDendro = c(0,4,1,2), marHeatmap = c(3,4,1,2), cex.lab = 0.8, yLabelsAngle = 90)
WGCNA-9

拆分树状图和热图图

# Plot the dendrogram
# sizeGrWindow(6,6);
par(cex = 1.0)
plotEigengeneNetworks(MET, "Eigengene dendrogram", marDendro = c(0,4,2,0),
plotHeatmaps = FALSE)
# Plot the heatmap matrix (note: this plot will overwrite the dendrogram plot)
par(cex = 1.0)
plotEigengeneNetworks(MET, "Eigengene adjacency heatmap", marHeatmap = c(3,4,2,2),plotDendrograms = FALSE, yLabelsAngle = 90)

版权声明：
作者：lichengxin
链接：https://www.techfm.club/p/55066.html
来源：TechFM
文章版权归作者所有，未经允许请勿转载。

THE END
分享
二维码
< <上一篇
下一篇>>