美文网首页组学学习
2023-12-02转录组分析——WGCNA

2023-12-02转录组分析——WGCNA

作者: 麦冬花儿 | 来源:发表于2023-12-01 08:56 被阅读0次

    Weighted Gene Co-Expression Network Analysis(加权基因共表达网络分析,简称WGCNA)可鉴定表达模式相似的基因集合(module),解析基因集合与样品表型之间的联系,绘制基因集合中基因之间的调控网络并鉴定关键调控基因。
    WGCNA的流程主要分为:输入数据——构建共表达网络——划分模块——模块与性状关联分析——模块之间的关联分析——模块中Hub基因的鉴定。此处输入的数据是经过其他方法筛出来的,比如之前TCseq分析中,对照和处理之间有差异的Cluster集群。具体原理网上有很多,这里就不多说,直接开搞。

    1、文件准备

    首先,就是准备两个文件:
    1、基因表达矩阵文件,即横轴是不同样本以及生物重复ID,纵轴是基因ID,之后是基因的FPKM值,在此将该文件命名为“FPKM.txt”。

    图片.png

    2、性状数据文件,横轴是不同性状的名称,纵轴是对应样品的名称,格式必须与上个文件保持一致,在此将该文件命名为“trait.txt”。


    图片.png

    2、构建相关性矩阵和邻接矩阵

    通过R语言中的WGCNA包开始分析,要注意的是WGCNA依赖的包较多,根据系统提示,逐一安装其他包。

    #安装WGCNA包
    install.packages('BiocManager') 
    BiocManager::install('WGCNA')
    library(WGCNA)  #若加载失败,安装提示中依赖的包
    #基因表达值矩阵
    gene <- read.delim('FPKM.txt', row.names = 1, check.names = FALSE)
    #只保留平均表达值在 1 以上的基因
    gene <- subset(gene, rowSums(gene)/ncol(gene) >= 1)
    #转置矩阵
    gene <- t(gene)
    

    接着需要我们确定一个软阈值powers(soft threshold),以建立邻接矩阵并根据连通度使基因分布符合无尺度网络。

    #power 值选择
    powers <- 1:20
    sft <- pickSoftThreshold(gene, powerVector = powers, verbose = 8)
     
    #拟合指数与 power 值散点图
    par(mfrow = c(1, 2))
    plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2], type = 'n', 
        xlab = 'Soft Threshold (power)', ylab = 'Scale Free Topology Model Fit,signed R^2', 
        main = paste('Scale independence'))
    text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
        labels = powers, col = 'red');
    abline(h = 0.80, col = 'red')  #R方的值可以自己确定,一般>0.8最好
    #平均连通性与 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, col = 'red')
    
    图片.png

    左图纵坐标是R平方的值,R方越接近1,表明该网络就越接近无尺度网络,通常要求R方的值>0.8,负值是因为乘slope列数值的负方向。软阈值一般选择R方大于0.8时的power值,即10。如果不确定,可以在代码中输入“sft$powerEstimate”让系统自身估计。


    图片.png

    可以看到系统评估的是13,emmm...当然也要综合上面右图——纵坐标是平均连通性,该值随β的增加而降低。那么接下来就用系统评估的最佳阈值13构建无标度网络。

    3、构建共表达网络

    #估计的最佳 power 值
    powers <- sft$powerEstimate
     
    #获得 TOM 矩阵
    adjacency <- adjacency(gene, power = powers)
    tom_sim <- TOMsimilarity(adjacency)
    rownames(tom_sim) <- rownames(adjacency)
    colnames(tom_sim) <- colnames(adjacency)
    tom_sim[1:6,1:6]
    #TOM 相异度 = 1 – TOM 相似度
    tom_dis  <- 1 - tom_sim
    #层次聚类树
    geneTree <- hclust(as.dist(tom_dis), method = 'average')
    plot(geneTree, xlab = '', sub = '', main = 'Gene Dendrogram',
        labels = FALSE, hang = 0.04)
    #使用动态剪切树挖掘模块
    minModuleSize <- 30  #模块基因数目
    dynamicMods <- cutreeDynamic(dendro = geneTree, distM = tom_dis,
        deepSplit = 2, pamRespectsDendro = FALSE, minClusterSize = minModuleSize)
    
    table(dynamicMods)   #查看共划分多少个模块
    
    图片.png

    该步骤是对基因进行聚类,每条线代表一个基因,相似的基因被聚到一个分支。输入“table(dynamicMods)”后可以看到,将这些基因划分为14个模块。


    图片.png

    接下来,为这些模块增添颜色,并通过颜色名称命名不同的模块。

    #模块颜色指代
    dynamicColors <- labels2colors(dynamicMods)
    table(dynamicColors)
     
    plotDendroAndColors(geneTree, dynamicColors, 'Dynamic Tree Cut',
        dendroLabels = FALSE, addGuide = TRUE, hang = 0.03, guideHang = 0.05,
        main = 'Gene dendrogram')
    
    图片.png

    4、共表达拓扑热图

    
    #基因表达聚类树和共表达拓扑热图
    plot_sim <- -(1-tom_sim)
    diag(plot_sim) <- NA
    TOMplot(plot_sim, geneTree, dynamicColors,
        main = 'Network heatmap plot')
    
    图片.png

    在该热图中,基因间表达相似度越高,颜色就会越深。

    #计算基因表达矩阵中模块的特征基因
    MEList <- moduleEigengenes(gene, colors = dynamicColors)
    MEs <- MEList$eigengenes
    head(MEs)[1:6]
    #表征模块间相似度
    ME_cor <- cor(MEs)
    ME_cor[1:6,1:6]
     
    #绘制聚类树
    METree <- hclust(as.dist(1-ME_cor), method = 'average')
    plot(METree, main = 'Clustering of module eigengenes', xlab = '', sub = '')
     
    #可以通过height 值确定一个合适的阈值作为剪切高度
    abline(h = 0.2, col = 'blue')
    abline(h = 0.25, col = 'red')
    #模块特征基因聚类树热图
    plotEigengeneNetworks(MEs, '', cex.lab = 0.8, xLabelsAngle= 90,
        marDendro = c(0, 4, 1, 2), marHeatmap = c(3, 4, 1, 2))
    #相关程度高于 0.75 的模块将合并到一起
    merge_module <- mergeCloseModules(gene, dynamicColors, cutHeight = 0.25, verbose = 3)
    mergedColors <- merge_module$colors
    table(mergedColors)
     #基因表达和模块聚类树
    plotDendroAndColors(geneTree, cbind(dynamicColors, mergedColors), c('Dynamic Tree Cut', 'Merged dynamic'),
        dendroLabels = FALSE, addGuide = TRUE, hang = 0.03, guideHang = 0.05)
    

    5、模块与性状关联

    接下来就是性状数据与基因模块之间的关联分析。

    #导入trait性状数据
    trait <- read.delim('trait.txt', row.names = 1, check.names = FALSE)
     
    #使用上一步新组合的共表达模块的结果
    module <- merge_module$newMEs
     
    #基因共表达模块和表型的相关性分析
    moduleTraitCor <- cor(module, trait, use = 'p')
    moduleTraitCor[1:6,1:6]  #相关矩阵,根据自己的性状个数进行修改
     
    #相关系数的 p 值矩阵
    moduleTraitPvalue <- corPvalueStudent(moduleTraitCor, nrow(module))
     
    #相关图绘制
    textMatrix <- paste(signif(moduleTraitCor, 2), '\n(', signif(moduleTraitPvalue, 1), ')', sep = '')
    dim(textMatrix) <- dim(moduleTraitCor)
     
    par(mar = c(5, 10, 3, 3))
    labeledHeatmap(Matrix = moduleTraitCor, main = paste('Module-trait relationships'),
        xLabels = names(trait), yLabels = names(module), ySymbols = names(module),
        colorLabels = FALSE, colors = blueWhiteRed(50), cex.text = 0.7, zlim = c(-1,1),
        textMatrix = textMatrix, setStdMargins = FALSE)
    
    图片.png
    每一个模块中,上面的数字代表皮尔逊相关系数,下面(括号)里的值代表显著性P值。该图表示哪些模块中的基因可能与处理表型密切相关。从图中不难发现,“green”模块“MTA10”时期存在较为显著的正相关,表示该模块中的基因参与了MTA处理引起的表型差异,进一步获取该模块中的基因(ps:该数据并不好,相关性系数不是太高,只有0.7,一般最好在0.9左右)。

    6、筛选关键基因集

    #基因与模块的对应关系列表
    gene_module <- data.frame(gene_name = colnames(gene), module = mergedColors, stringsAsFactors = FALSE)
    head(gene_module)
    #“green”模块内的基因名称
    gene_module_select <- subset(gene_module, module == 'green')$gene_name
     
    #“green”模块内基因在各样本中的表达值矩阵
    gene_select <- t(gene[,gene_module_select])
     
    #“green”模块内基因的共表达相似度
    tom_select <- tom_sim[gene_module_select,gene_module_select]
    #选择 green 模块内的基因
    gene_green <- gene[ ,gene_module_select]
     
    #基因的模块成员度计算
    geneModuleMembership <- signedKME(gene_green, module['MEgreen'], outputColumnName = 'MM')
    MMPvalue <- as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nrow(module)))
     
    #各基因表达值与表型的相关性分析
    geneTraitSignificance <- as.data.frame(cor(gene_green, trait['MTA10'], use = 'p'))
    GSPvalue <- as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nrow(trait)))
     
    #选择显著(p<0.01)、高 green 模块成员度(MM>=0.8),与 MTA10表型高度相关(r>=0.8)的基因
    geneModuleMembership[geneModuleMembership<0.8 | MMPvalue>0.01] <- 0
    geneTraitSignificance[geneTraitSignificance<0.8 | MMPvalue>0.01] <- 0
     
    select <- cbind(geneModuleMembership, geneTraitSignificance)
    select <- subset(select, geneModuleMembership>=0.8 & geneTraitSignificance>=0.8)
    head(select)
    nrow(select)#查看相关系高的基因个数
    
    图片.png

    最终筛选到了16个候选基因(右侧“Environment”中点击“select”即可查看),进一步结合注释信息确认其是否是自己想要的hub基因。WGCNA分析过程较为繁琐,一环扣一环,要确保每一步正确,再进行下一步操作。

    相关文章

      网友评论

        本文标题:2023-12-02转录组分析——WGCNA

        本文链接:https://www.haomeiwen.com/subject/gwlcgdtx.html