WGCNA分析（4-21）

于 2025-05-03 19:44:20 发布
阅读量716 收藏 19
点赞数 14
CC 4.0 BY-SA版权
文章标签： 机器学习 人工智能
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/weixin_42630015/article/details/144703956

 R包加载

#安装R包
library(WGCNA)
library(ggplot2)
library(reshape2)
library(stringr)
library(dplyr)

读入基因表达数据

#1.基因表达数据整理

#1.1数据读入
PAAD_FPKM_data<-fread("C:/Users/hp/Desktop/t/TCGA-PAAD.star_fpkm.tsv.gz",header = T,sep = "\t",data.table = F)

#1.2基因ID转换
PAAD_gene<-fread("C:/Users/hp/Desktop/t/gencode.v36.annotation.gtf.gene.probemap",header = T,sep = "\t",data.table = F)
PAAD_gene<-PAAD_gene[,c(1,2)]
PAAD_FPKM_data_gene<-merge(PAAD_gene,PAAD_FPKM_data,by.x="id",by.y="Ensembl_ID")

PAAD_FPKM_data数据框

PAAD_gene数据框

数据预处理

去除重复基因

dim(PAAD_FPKM_data_gene)

#1.3去除重复基因
sum(duplicated(PAAD_FPKM_data_gene$gene))

#distinct用于去除数据框中的重复行。
PAAD_FPKM_data_gene<-distinct(PAAD_FPKM_data_gene,gene,.keep_all=T)    #gene为指定了用于识别重复行的列名；即使某些行在gene列上是重复的，keep_all=T也会保留这些行的所有列信息
dim(PAAD_FPKM_data_gene)

 行列转制

PAAD_FPKM_data_gene<-column_to_rownames(PAAD_FPKM_data_gene,"gene")
PAAD_FPKM_data_gene<-PAAD_FPKM_data_gene[,-1]
#1.4数据行列转换
PAAD_FPKM_data_gene <- data.frame(t(PAAD_FPKM_data_gene))
dim(PAAD_FPKM_data_gene)
PAAD_FPKM_data_gene[1:3,1:3]

剔除正常样本，保留肿瘤样本

table(str_sub(rownames(PAAD_FPKM_data_gene),14,15))

group_list= ifelse(as.numeric(str_sub(rownames(PAAD_FPKM_data_gene),14,15)) < 10,'tumor','normal')
group_list= factor(group_list,levels = c("normal","tumor"))
table(group_list)
gene_data <- PAAD_FPKM_data_gene %>% filter(group_list == "tumor")

剔除缺失值、权重低于阈值的条目以及零方差基因 

#1.5.2剔除缺失值、权重低于阈值的条目以及零方差基因
gsg = goodSamplesGenes(gene_data, verbose = 3)   ##用于检查数据集中的缺失值、权重低于阈值的条目以及零方差基因
if (!gsg$allOK) {
  # 输出需要删除的基因和样本名称
  if (sum(!gsg$goodGenes)>0)
    printFlush(paste("Removing genes:", paste(names(gene_data)[!gsg$goodGenes], collapse = ", ")))
  if (sum(!gsg$goodSamples)>0)
    printFlush(paste("Removing samples:", paste(rownames(gene_data)[!gsg$goodSamples], collapse = ", ")))
  # 删除
  gene_data = gene_data[gsg$goodSamples, gsg$goodGenes]
}
dim(gene_data)

剔除离群数据 

#1.5.3剔除离群数据

# 方法1.查看是否有离群样品，hclusts聚类算法, dist计算基因之间的距离
sampleTree = hclust(dist(gene_data), method = "average")

par(cex = 0.2)  # 控制图片中文字和点的大小
par(mar =c(1,5,2,1))    # 设置图形的边界，下，左，上，右的页边距
plot(sampleTree, main ="Sample clustering to detectoutliers",sub="", xlab="", cex.lab = 1.5,cex.axis= 1.5, cex.main = 2)

# 可以手动删除软阈值，也可以使用以下办法
# 绘制阈值切割线
abline(h = 100,col="red") # 高度100，颜色红色
# 确定阈值线下的集群
clust = cutreeStatic(sampleTree, cutHeight = 100, minSize = 10)
# 以高度100切割，要求留下的最少为10个
table(clust)    # 查看切割后形成的集合
# clust1包含想要留下的样本.
keepSamples = (clust==1)    # 将clust序号为1的放入keepSamples
dataExpr = dataExpr[keepSamples, ]      
# 将树中内容放入矩阵dataExpr中，因为树中剩余矩阵不能直接作为矩阵处理
nGenes = ncol(dataExpr)
nSamples = nrow(dataExpr)
dim(dataExpr)

sampleTree 数据结构：

merge:描述每一步聚类时两个子簇（或观测点）的合并过程。 

  • 负值表示原始观测点（如 -3 表示第 3 个样本）。

  • 正值表示之前合并生成的新簇（如 2 表示第 2 步生成的簇）。

  • merge:
     [,1] [,2]
[1,]   -1   -2   # 第1步：合并样本1和2，生成新簇1
[2,]   -3    1   # 第2步：合并样本3和新簇1，生成新簇2
[3,]   -4   -5   # 第3步：合并样本4和5，生成新簇3

height:记录每次合并时两个子簇之间的距离值。（以merge合并顺序）
height:
[1] 0.5 1.2 2.0  # 第1步合并距离0.5，第2步1.2，第3步2.0

order

  • 含义: 树状图中叶子节点（观测点）的排列顺序，用于绘制无交叉的树状图。

  • order:
[1] 3 1 2 4 5     # 树状图从左到右显示的顺序为：观测点3 → 1 → 2 → 4 → 5

    此数据为例

如上数据表示第10个样本与第134个样本先合并，然后第139个样本与第1个新簇合并，即第10与第134合并后再与第139合并。

此树从左到右分别是第118、13。。。

第10：TCGA-F2-7276；第134：TCGA-F2-7273；第139：TCGA-HZ-8002

# 方法2.此外，可以通过计算连接矩阵的方式来评估是否有离群值
A = adjacency(t(gene_data),type = "distance")
k = as.numeric(apply(A,2,sum)) - 1
Z.k = scale(k)
# 设置离群值的阈值
thresholdZ.k = -5 #often -2.5
# 离群值为红色，否则为黑色
outlierColor = ifelse(Z.k < thresholdZ.k, "red","black")

sampleTree = hclust(as.dist(1-A), method = "average")
plotDendroAndColors(sampleTree,groupLabels = names(outlierColor),
                    colors = outlierColor,main = "Sample dendrogram and trait heatmap")
                    
# 方法2.end

确定软阈值 

powers = c(c(1:10), seq(from = 12, to=30, by=2))
# 调用网络拓扑分析函数
#sft = pickSoftThreshold(PAAD_FPKM_data_gene, powerVector = powers, networkType = type, verbose = 5)

sft <- pickSoftThreshold(as.matrix(gene_data),
                         powerVector = powers, 
                         networkType = "signed")

 # 绘图
sizeGrWindow(9, 5)
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 = "Scale independence")
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     labels=powers,cex=0.9,col="red")
# 筛选标准：R square=0.85
abline(h=0.85,col="red")
# 软阈值与平均连通性
plot(sft$fitIndices[,1], sft$fitIndices[,5],
     xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
     main = "Mean connectivity")
text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=0.9,col="red")

power = sft$powerEstimate

sft数据结构：

 以下代码含义

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 = "Scale independence") 

  • -sign(sft$fitIndices[,3])*sft$fitIndices[,2]:

    • sft$fitIndices[,3] 是斜率（slope列）；sign() 用于判断数值的正负性，返回 1（正数）、0（零）或 -1（负数）。

    • -sign(slope) 的目的是将负斜率转换为正值（理想无标度网络的斜率为负，因此此处实际取正值）。

    • 乘以 SFT.R.sq（第2列）得到最终的拟合指数。

  • type = "n": 不绘制任何点或线，仅初始化画布。

构建共表达网络 

nGenes <- ncol(gene_data)

net = blockwiseModules(gene_data, power = power, maxModuleSize = nGenes,
                       TOMType = "unsigned", minModuleSize = 25, 
                       reassignThreshold = 0, 
                       mergeCutHeight = 0.2, numericLabels = TRUE,
                       pamRespectsDendro = FALSE, saveTOMs = TRUE,
                       verbose = 3,randomSeed = 1234)
 # 这一步需要运算TOM值，运行会比较慢，并且需要确保内存足够。
table(net$colors)

net文件

①colors:每个基因所属模块的标签(126个模块)

②unmergedColors为未合并前基因所属模块

③ME：每个模块的特征基因（eigengene），这些特征基因的表达（179行 126列）

④dendrograms：基因的层次聚类树

⑤blockGenes

 计算MM与GS

  1. MM (Module Membership)：

    • 模块成员关系。它是一个量化度量，表示基因表达谱与模块特征基因（eigengene）之间的相关性。MM值本质上是一个相关系数，用来衡量一个基因与特定模块的关联程度。如果一个基因的MM值接近0，说明该基因与模块关系不大；如果MM值的绝对值接近1，说明该基因与模块高度相关。

  2. GS (Gene Significance)：

    • 基因显著性。它量化了单个基因与感兴趣的表型特征（如体重、疾病状态等）之间的关联程度，定义为基因表达数据与表型数据之间的相关性的绝对值。GS值越大，表示基因表达与表型特征的相关性越强。
# 1.MM计算
##提取模块名称
modNames = substring(names(MEs), 3)  # names与colname指列名

#计算基因与模块的相关系数矩阵（计算皮尔斯系数）
geneModuleMembership = as.data.frame(cor(gene_data, MEs,
                                           use = "p",method = "spearman")) ##共54374行（基因），126列（模块） 

MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples)) （dim数据同上）

 

#修改列名
names(geneModuleMembership) = paste("MM", modNames, sep="")  # names与colname指列名
names(MMPvalue) = paste("p.MM", modNames, sep="")  # 修改列名


#2.计算GS数据
geneTraitSignificance <- as.data.frame(cor(gene_data,datTraits$ajcc_pathologic_stage.diagnoses,use = "p"))  # 计算基因与表型的相关系数矩阵（共54374行（基因），1列（临床特征） ）
GSPvalue <- as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance),nSamples))  # 计算相关系数矩阵的p值（共54374行（基因），1列（临床特征） ）
names(geneTraitSignificance)<- paste("GS." , "ajcc_pathologic_stage.diagnoses",sep = "")  # 修改列名
names(GSPvalue)<-paste("GS." , "ajcc_pathologic_stage.diagnoses",sep = "")  # 修改列名

geneTraitSignificance数据框

 

GSPvalue数据框

后两行重命名是将“GS.ajcc_pathologic_stage.diagnoses”设置为数据列名。

#画MM-GS图

selectModule<-c("brown","red")  # 选择要绘制的模块（这里只选择了一个）
#selectModule <- modNames  # 批量作图

  
par(mfrow=c(ceiling(length(selectModule)/2),2))  # 设置绘图区域，批量作图开始
fs <- list() 
for(module in selectModule){
    column <- match(module,modNames)  # 找到当前模块在geneModuleMembership中的列号
    print(module)  # 输出当前处理的模块名称
    moduleGenes <- moduleColors==module  # 找到属于当前模块的基因

    # 绘制散点图，横轴为基因在当前模块中的模块内连通性，纵轴为基因与表型的相关系数
verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
                       abs(geneTraitSignificance[moduleGenes, 1]),
                       xlab = paste("Module Membership in", module, "module"),
                       ylab = paste("Gene significance for", module, "module"),
                       main = paste("Module membership vs. gene significance\n"),
                       cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)

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

rownames(f) = names(moduleLabels[moduleGenes])

fs[[module]] <- f
  }

 ***将临床性状数据由“字符”改为“数字”

#读入临床数据
PAAD_sur<-fread("C:/Users/hp/Desktop/t/TCGA-PAAD.survival.tsv.gz",header = T,sep = "\t",data.table = F)
PAAD_clin<-fread("C:/Users/hp/Desktop/t/TCGA-PAAD.clinical.tsv.gz",header = T,sep = "\t",data.table = F)
PAAD_clin<-PAAD_clin[,c(1,3,8,11,13,14,45,48)]

#将字符型变量替换为字符型变量
PAAD_clin$disease_type <- ifelse(PAAD_clin$disease_type == "Adenomas and Adenocarcinomas", 1,
                          ifelse(PAAD_clin$disease_type == "Cystic, Mucinous and Serous Neoplasms", 2,
                          ifelse(PAAD_clin$disease_type == "Ductal and Lobular Neoplasms", 3, 
                          ifelse(PAAD_clin$disease_type == "Epithelial Neoplasms, NOS", 4,                                          5))))

PAAD_clin$alcohol_history.exposures <- ifelse(PAAD_clin$alcohol_history.exposures == "Not Reported", 0,
                                       ifelse(PAAD_clin$alcohol_history.exposures == "No", 1,
                                       ifelse(PAAD_clin$alcohol_history.exposures == "Yes", 2, 
                                       3)))

PAAD_clin$gender.demographic <- ifelse(PAAD_clin$gender.demographic == "female", 0,
                                ifelse(PAAD_clin$gender.demographic == "male", 1,
                                2))

PAAD_clin$vital_status.demographic <- ifelse(PAAD_clin$vital_status.demographic == "Alive", 0,
                                      ifelse(PAAD_clin$vital_status.demographic == "Dead", 1,
                                      2))

PAAD_clin$ajcc_pathologic_stage.diagnoses <- ifelse(PAAD_clin$ajcc_pathologic_stage.diagnoses == "Stage I", 1,
                                      ifelse(PAAD_clin$ajcc_pathologic_stage.diagnoses == "Stage IA", 1,
                                      ifelse(PAAD_clin$ajcc_pathologic_stage.diagnoses == "Stage IB", 2,                                        ifelse(PAAD_clin$ajcc_pathologic_stage.diagnoses == "Stage IIA", 3,
                                      ifelse(PAAD_clin$ajcc_pathologic_stage.diagnoses == "Stage IIB", 4,                                       ifelse(PAAD_clin$ajcc_pathologic_stage.diagnoses == "Stage III", 5,                                       ifelse(PAAD_clin$ajcc_pathologic_stage.diagnoses == "Stage IV", 6,                                        0)))))))

PAAD_clin$tissue_or_organ_of_origin.diagnoses <- ifelse(PAAD_clin$tissue_or_organ_of_origin.diagnoses == "Head of pancreas", 1,
                                      ifelse(PAAD_clin$tissue_or_organ_of_origin.diagnoses == "Body of pancreas", 2,
                                      ifelse(PAAD_clin$tissue_or_organ_of_origin.diagnoses == "Tail of pancreas", 3,
                                      ifelse(PAAD_clin$tissue_or_organ_of_origin.diagnoses == "Overlapping lesion of pancreas", 4,   
                                                 0))))

***将临床数据的样本名排列顺序调整为与基因表达数据（gene_data）样本名排列相同

agenes<-setdiff(PAAD_clin[["sample"]],rownames(gene_data))
#在离群值筛选时可能会删除部分样品，在这里需要将对应的表型数据去除

##allTraits=data.frame(PAAD_clin[-which(PAAD_clin[["sample"]]==agenes),])
allTraits = data.frame(PAAD_clin[-which(PAAD_clin[["sample"]] %in% agenes), ])
matchSamples=rownames(gene_data)
traitRows=match(matchSamples,allTraits$sample)
datTraits=allTraits[traitRows,-1]
rownames(datTraits)=allTraits[traitRows,1]

 1.agenes<-setdiff(PAAD_clin[["sample"]],rownames(gene_data))

  • setdiff(x, y) 是 R 中的一个函数，用于计算集合 x 中有但集合 y 中没有的元素。

  • 在这段代码中，setdiff(PAAD_clin[["sample"]], rownames(gene_data)) 的作用是找出 PAAD_clin[["sample"]] 中有但 rownames(gene_data) 中没有的元素。

2.allTraits = data.frame(PAAD_clin[-which(PAAD_clin[["sample"]] %in% agenes), ])

  1. which(PAAD_clin[["sample"]] == agenes):

    • which() 函数返回逻辑向量中为 TRUE 的位置索引。

    • 也就是说，它会返回 PAAD_clin[["sample"]] 中与 agenes 匹配的行号。

  2. -which(PAAD_clin[["sample"]] == agenes):

    • 负号 - 用于从索引中排除这些行。

196-17=179，所以195个样本是有问题的

3.traitRows=match(matchSamples,allTraits$sample)

示例：

allTraits <- data.frame(
             sample = c("sample1", "sample2", "sample3", "sample4"),
              age = c(45, 50, 55, 60),
              status = c("A", "B", "A", "B")
)

matchSamples <- c("sample2", "sample4", "sample5")

traitRows = match(matchSamples, allTraits$sample)

traitRows

[1]  2  4 NA

***画聚类图（样本与临床信息） 

#画聚类图
sampleTree2=hclust(dist(gene_data),method="average")
#画表型的热图
#将特征转换为颜色表示：白色表示低，红色表示高，灰色表示缺少条目
traitColors=numbers2colors(datTraits,signed=FALSE)
#绘制出树状图和热图,为样本聚类与表型的相关性
plotDendroAndColors(sampleTree2,traitColors,groupLabels=names(datTraits),main="Sample dendrogram and trait heatmap")

 1.dist(gene_data)结果：计算 gene_data 中每一行之间的距离，并返回一个距离对象

2.sampleTree2=hclust(dist(gene_data),method="average")结果：hclust() 函数对距离矩阵进行层次聚类分析

3.traitColors=numbers2colors(datTraits,signed=FALSE)：将临床数据数值改为颜色

MEs_colpheno = orderMEs(cbind(MEs_col,datTraits))
par(cex = 0.5)
plotEigengeneNetworks(MEs_colpheno,
                      "Eigengene adjacrncy heatmap",
                      marDendro = c(3,3,2,4),marHeatmap = c(3,4,2,2),
                      plotDendrograms = T,
                      xLabelsAngle = 90)

              MEs_col数据框                                               datTraits数据框

1.cbind（）是以行名为依据，将两个数据合并，

2.orderMEs（）将给定的（特征）向量重新排序，使得相似的向量（依据相关性（correlation））彼此相邻。

                 MEs_colpheno数据框（仅列排列顺序改变，其他不变，包括数值）

3.plotEigengeneNetworks函数会计算模块与模块的相关性，并绘制成热图。

# 用彩色标签重新计算MEs,在给定的单个数据集中计算模块的模块本征基因
MEs0 = moduleEigengenes(gene_data, moduleColors)$eigengenes
# 对给定的（特征）向量进行重新排序，以使相似的向量（通过相关性度量）彼此相邻
MEs = orderMEs(MEs0)
# 计算module的ME值与表型的相关系数
moduleTraitCor = cor(MEs, datTraits, use = "p")
nSamples = nrow(gene_data)
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)

1.moduleColors数据框

2.gene_data数据框

3.moduleEigengenes(gene_data, moduleColors)结果

4.moduleEigengenes(gene_data, moduleColors)$eigengenes结果

5.MEs = orderMEs(MEs0)结果：（列的排列顺序改变，余不变）

6.moduleTraitCor = cor(MEs, datTraits, use = "p")结果：计算各模块与临床性状的相关系数

7.moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)计算6的相关系数的P值

8、

在WGCNA分析流程中，blockwiseModules函数返回的MEs（Module Eigengenes）与后续调用moduleEigengenes生成的MEs0可能在表面上结果一致，但二者的设计目的和使用场景存在差异。以下是关键原因：

1. 计算阶段的差异

  • blockwiseModules中的MEs：
    该函数在模块识别过程中自动计算模块特征基因（MEs）。这一过程基于算法内部定义的模块划分结果（即moduleColors），且可能涉及模块合并、动态剪切等步骤。返回的MEs反映的是最终优化后的模块结构对应的特征基因。

  • moduleEigengenes的MEs0：
    该函数是独立工具，允许用户在任意阶段基于当前模块标签（如手动调整后的moduleColors）重新计算MEs。例如，若后续合并模块或调整颜色后，需用此函数更新MEs。

2. 动态调整模块的需求

  • 若在初始模块识别后，用户通过mergeCloseModules等函数合并模块，或手动修改moduleColors，原始net$MEs将不再与新的模块结构匹配。此时必须调用moduleEigengenes，基于更新后的模块标签重新计算MEs，确保后续分析（如模块-性状关联）的准确性。

3. 数据一致性的验证

  • 当MEs与MEs0结果一致时，表明模块划分未发生后续变动，且数据预处理一致。这种验证有助于确认模块识别的稳定性，避免因参数设置或数据输入错误导致的偏差。

9、相关系数大小的意义

# 相关性及其p值可视化
textMatrix = paste(signif(moduleTraitCor, 2), "\n(",signif(moduleTraitPvalue, 1), ")", sep = "")
dim(textMatrix) = dim(moduleTraitCor)
par(mfrow = c(1,1), cex = 1.0, mar = c(6, 8, 2, 1))
labeledHeatmap(Matrix = moduleTraitCor,xLabels = names(datTraits),yLabels = names(MEs),ySymbols = names(MEs),
               colorLabels = FALSE,colors = greenWhiteRed(50),textMatrix = textMatrix,setStdMargins = FALSE,
               cex.text = 0.4,zlim = c(-1,1),
               main = paste("Module-trait relationships"))

gene_data表格：

 

MEs表格：

geneModuleMembership表格如下：

MMPvalue表格如下：

geneTraitSignificance表格如下：

MM：
定义

  • MM（模块成员资格）衡量某个基因与所属模块的 相关性强度，即该基因在模块中的表达模式与模块整体表达模式的一致性。

  • 取值范围：[-1, 1]，绝对值越接近1，基因与模块的关联越强。

应用

  • 筛选模块核心基因：高|MM|值的基因（如|MM| > 0.8）可能是模块的枢纽基因（hub genes）。

 GS：

定义

  • GS（基因显著性）衡量某个基因与 外部表型特征（如疾病状态、生存时间、实验处理等）的 相关性强度。

取值范围

  • 若为相关系数：[-1, 1]，绝对值越接近1，基因与表型关联越强。

应用

  • 识别表型相关基因：高|GS|或低p值的基因可能与表型直接相关。

  • 模块与表型关联分析：结合模块的ME与表型的相关性（Module-Trait Analysis），筛选关键模块。

藍色飛揚
关注
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值