一步完成WGCNA

2023-12-30 11:32
文章标签 完成 一步 wgcna

本文主要是介绍一步完成WGCNA，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

OneStepWGCNA
今天发现TBtools上有了个神器 Rserver 。这意味我们可以把别人写好的R脚本直接拿过来运行。于是用OneStepWGCNA这个插件试验一下。本文主要记录试验过程和可能遇到的一些问题，主要参考『TBtools-plugin』OneStepWGCNA插件这篇文章。想要了解WGCNA的具体含义和流程可以参考：
WGCNA分析，简单全面的最新教程
一文学会WGCNA分析
WGCNA官方网站WGCNA算法研究笔记
OneStepWGCNA 的准备文件和详细参数

虽然题目有点标题党，但这个插件真的方便，真的感谢这个插件的作者。

工具准备

TBtools
Rserver: 需要在TBtools的用户群里下载，一共181M，相当于一个独立的R环境，具体安装方法参考Rserver 插件 for TBtools。用户群可以在公众号：生信药丸里找到。
OneStepWGCNA插件：可以在TBtools的plugin store下载到。建议先下载plugin store里的plugin store at high speed。因为在plugin store里下载插件太慢了，太慢了，太慢了。
#数据准备
OneStepWGCNA 的输入文件是基因表达矩阵和表型数据（或样品分类数据）。这里用了RED数据库里的284个日本晴的不同部位，不同处理的表达矩阵(FPKM数据)。样品分类数据用每个样品取的组织。
表达矩阵
样品分类数据

第一次跑的时候遇到的坑

报错信息
第一次包发现GenomeInfoDb这个包没下下来，R文件里也没有申明需要这个包。我在TBtools的插件目录下OneStepWGCNA的R脚本里加上了：

if (!require(‘GenomeInfoDb’)) install.packages(‘GenomeInfoDb’)

R文件
再跑一次就成功了
这个TBtools插件文件夹位置一般在C:\Users\names\ .TBtools.Plugin\OneStepWGCNA

OneStepWGCNA 参数详解

readcount or fpkm 是表达矩阵，以Tab为分割符，行是样本，列是基因。
Trait data 是表型文件或分类文件，可以有多个表型，可以是连续性也可以是离散。如果要把多分类性状和连续性状都混到一起，可以把多分类性状转为one-hot编码，这个用下excel就能完成。
one-hot编码
Expression data 是表达矩阵数据类型，count数据是不能直接用来做WGCNA的
Normalization method 是表达矩阵转换的方式，如果是count数据可以选varianceStabilizingTransformation或log10(CPM+1)，如果是FPKM可以选rawFPKM，不经过转换，和取对数。建议都试试，看看哪个效果更好。
RcCutoff 是用来除噪声的，我看原作者的教程里大概意思是count推荐10，FPKM推荐0.2，当然，这个参数还是要随着测序深度和具体生物学问题而改变的。
samplePerc 和上面的参数配对的，在0-1之间，如果是0.5，RcCutoff 是0.2 代表想要 50%的样品的FPKM>0.2。
RemainGeneNum 是保留多少基因进行WGCNA，写0或小于0就是啥都没了，会报错。
mergeCutHeight 是合并模块的阈值， 0.25的意思就是把相关系数大于0.75的模块合并。
minModuleSize 是指最少一个模块要包含30个基因。

输出结果

Sample_clustering 样品聚类结果
样品聚类结果
这里如果出现明显的离群样本，记住样本的名字，在表达矩阵和trait矩阵中删掉对应样本再跑一次
软阈值筛选
绿线为0.8，这里没有筛的很好的阈值，就选了经验阈值14来做power进行后续计算。
模块的邻接矩阵特征向量相关性
一共找到了10个模块，该结果主要反应模块间的相关性。
模块划分
表型与模块的关联
发现棕色模块和根的相关性最高，其他种子，花药，穗的部位都没有找到相关性高的模块。一个可能是由于这几个部位的样本量少，另一方面本身不同处理对表达量的影响也很大。这里只是试一试，并不能反应真实的科学结论。
基因对应的模块
可以拿表型关联上的模块做后续的分析，如GO，KEGG啥的，巧的是富集分析也能用TBtools做。

最后再唠叨两句

我发现最后结果里还有很多信息没有被存下来，我想看选几个基因画TOM图，或者画个网络图，或者是图的字体太小了，看不清想重新画，该咋办？
最粗暴的办法还是把整个网络保存下来，想要啥，自己整。
修改了下blockwiseModules的参数，把TOM和net都保存下来
以下代码主要来源于 WGCNA分析，简单全面的最新教程

library(WGCNA)
#把存入的Rdata都读进来，总能用上的 \xk one_traits是我设的Title name,要改
load("00.one_traits.datatraitbase.Rdata")
load("00.one_traits.net.Rdata")
load("00.one_traits.datatraitbase.Rdata")
load("00one_traits.Modular.Rdata")
load("one_traits.tom-block.1.RData")
Title<-"one_traits"
###TOM图
moduleColors = labels2colors(net$colors)
TOM <- as.matrix(TOM)
dissTOM = 1-TOM
# Transform dissTOM with a power to make moderately strong 
# connections more visible in the heatmap
plotTOM = dissTOM^7
# Set diagonal to NA for a nicer plot
diag(plotTOM) = NA
# Call the plot function
# 这一部分特别耗时，行列同时做层级聚类
TOMplot(plotTOM, net$dendrograms, moduleColors, main = "Network heatmap plot, all genes")###导出Cytoscape 图格式
probes = colnames(datExpr)
dimnames(TOM) <- list(probes, probes)
# Export the network into edge and node list files Cytoscape can read
# threshold 默认为0.5, 可以根据自己的需要调整，也可以都导出后在
# cytoscape中再调整
cyt = exportNetworkToCytoscape(TOM,edgeFile = paste(Title, ".edges.txt", sep=""),nodeFile = paste(Title, ".nodes.txt", sep=""),weighted = TRUE, threshold = 0.5,nodeNames = probes, nodeAttr = moduleColors)## 模块内基因与表型数据关联
#关联样品性状的二元变量时，设置
robustY = ifelse(corType=="pearson",T,F)
if (corType=="pearsoon") {geneModuleMembership = as.data.frame(cor(datExpr, MEs_col, use = "p"))MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples))
} else {geneModuleMembershipA = bicorAndPvalue(datExpr, MEs_col, robustY=robustY)geneModuleMembership = geneModuleMembershipA$bicorMMPvalue   = geneModuleMembershipA$p
}
# 计算性状与基因的相关性矩阵
## 只有连续型性状才能进行计算，如果是离散变量，在构建样品表时就转为0-1矩阵。
traitData<-read.table(traitsfile,header=T,sep='\t') #traitsfile 表型文件
if (corType=="pearsoon") {geneTraitCor = as.data.frame(cor(datExpr, traitData, use = "p"))geneTraitP = as.data.frame(corPvalueStudent(as.matrix(geneTraitCor), nSamples))
} else {geneTraitCorA = bicorAndPvalue(datExpr, traitData, robustY=robustY)geneTraitCor = as.data.frame(geneTraitCorA$bicor)geneTraitP   = as.data.frame(geneTraitCorA$p)
}## Warning in bicor(x, y, use = use, ...): bicor: zero MAD in variable 'y'.
## Pearson correlation was used for individual columns with zero (or missing)
## MAD.# 最后把两个相关性矩阵联合起来,指定感兴趣模块进行分析
module = "brown"#感兴趣的模块
pheno = "Root" #感兴趣的性状
modNames = substring(colnames(MEs_col), 3)
# 获取关注的列
module_column = match(module, modNames)
pheno_column = match(pheno,colnames(traitData))
# 获取模块内的基因
moduleGenes = moduleColors == module# 与性状高度相关的基因，也是与性状相关的模型的关键基因
verboseScatterplot(abs(geneModuleMembership[moduleGenes, module_column]),abs(geneTraitCor[moduleGenes, pheno_column]),xlab = paste("Module Membership in", module, "module"),ylab = paste("Gene significance for", pheno),main = paste("Module membership vs. gene significance\n"),cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)

10000个基因的TOM
棕色模块与根的相关性

这篇关于一步完成WGCNA的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！


http://www.chinasem.cn/article/552691。

相关文章

Java 后端接口入参 - 联合前端VUE 使用AES完成入参出参加密解密

Java 后端接口入参 - 联合前端VUE 使用AES完成入参出参加密解密

加密效果： 解密后的数据就是正常数据： 后端：使用的是spring-cloud框架，在gateway模块进行操作 <dependency><groupId>com.google.guava</groupId><artifactId>guava</artifactId><version>30.0-jre</version></dependency> 编写一个AES加密
阅读更多...
一步一步将PlantUML类图导出为自定义格式的XMI文件

一步一步将PlantUML类图导出为自定义格式的XMI文件

一步一步将PlantUML类图导出为自定义格式的XMI文件 说明： 首次发表日期：2024-09-08PlantUML官网： https://plantuml.com/zh/PlantUML命令行文档： https://plantuml.com/zh/command-line#6a26f548831e6a8cPlantUML XMI文档: https://plantuml.com/zh/xmi
阅读更多...
2024年高教社杯数学建模国赛最后一步——结果检验-事关最终奖项

2024年高教社杯数学建模国赛最后一步——结果检验-事关最终奖项

2024年国赛已经来到了最后一天，有必要去给大家讲解一下，我们不需要过多的去关注模型的结果，因为模型的结果的分值设定项最多不到20分。但是如果大家真的非常关注的话，那有必要给大家讲解一下论文结果相关的问题。很多的论文，上至国赛优秀论文下至不获奖的论文并不是所有的论文都可以进行完整的复现求解，大部分数模论文都为存在一个灰色地带。         白色地带即认为所有的代码均可运行、公开
阅读更多...
如何完成本科毕业论文设计

如何完成本科毕业论文设计

完成本科毕业论文设计是一个系统性的工程，需要经过多个阶段的规划、执行和总结。以下是一个详细的步骤指南，帮助你顺利完成本科毕业论文设计。 ### 1. 选题与开题 - **选题**：选择一个有研究价值且你感兴趣的题目。与导师讨论，确保题目具有可行性和创新性。 - **开题报告**：撰写开题报告，包括研究背景、研究目的、研究内容、研究方法、预期成果等。 ### 2. 文献综述 - **文献检索**
阅读更多...
LabVIEW环境中等待FPGA模块初始化完成

LabVIEW环境中等待FPGA模块初始化完成

这个程序使用的是LabVIEW环境中的FPGA模块和I/O模块初始化功能，主要实现等待FAM（Field-Programmable Gate Array Module，FPGA模块）的初始化完成，并处理初始化过程中的错误。让我们逐步分析各部分的功能： 1. Wait for FAM Initialization框架 此程序框架用于等待I/O模块成功初始化。如果在5秒钟内模块没有完成配
阅读更多...
完成一个项目的流程

完成一个项目的流程

我自己总结的，有什么问题，请大家指点啊！ 1. 制定项目的周期。工具：project 2. 确定需求，设计界面。工具：Axure 3. 写需求文档。 4. 写接口文档。 5. 设计项目架构。工具：Visio 6. 做图。工具：ps 7. 编码。 8. 写测试用例。 9.测试。
阅读更多...
idea中配置Translation插件完成翻译功能

idea中配置Translation插件完成翻译功能

文章目录 idea下载插件配置有道云阿里云百度翻译开放平台 idea下载插件 idea中安装Translation插件 使用方法：右下角选择翻译引擎，鼠标选中想翻译的部分，右键翻译即可 之前一直用的微软的翻译，不需要配置，但是最近微软服务器总是抽风，无法使用，故打算配置一下国内的翻译服务。 配置 有道云 只有初始的一点额度，用完就要收费了，不推荐
阅读更多...
LiveQing视频点播流媒体RTMP推流服务功能-支持大疆等无人机RTMP推流支持OBS推流一步一步搭建RTMP视频流媒体服务示例

LiveQing视频点播流媒体RTMP推流服务功能-支持大疆等无人机RTMP推流支持OBS推流一步一步搭建RTMP视频流媒体服务示例

LiveQing支持大疆等无人机RTMP推流支持OBS推流一步一步搭建RTMP视频流媒体服务示例 1、流媒体服务搭建2、推流工具准备3、创建鉴权直播间4、获取推流地址5、配置OBS推流6、推流及播放7、获取播放地址7.1 页面查看视频源地址7.2 接口查询 8、相关问题8.1、大疆无人机推流花屏 9、RTMP推流视频直播和点播流媒体服务 1、流媒体服务搭建 Windows/Lin
阅读更多...
解决解压缩时的错误提示 “无法成功完成操作, 因为文件包含病毒或者潜在垃圾文件“

解决解压缩时的错误提示 “无法成功完成操作, 因为文件包含病毒或者潜在垃圾文件“

近期, 有一些朋友反馈在解压zip压缩包, 或者在安装软件的过程中出现了下面的错误提示: "无法成功完成操作, 因为文件包含病毒或者潜在垃圾文件" "Operation did not complete successfully because the file contains a virus or potentially unwanted software" 上述错误一般
阅读更多...
亚信安慧AntDB数据库与华为DPA数据保护一体机完成兼容性互认证，共筑数据安全与效率新高地

亚信安慧AntDB数据库与华为DPA数据保护一体机完成兼容性互认证，共筑数据安全与效率新高地

近日，湖南亚信安慧科技有限公司（简称“亚信安慧”）与华为技术有限公司（简称“华为”）完成了亚信安慧AntDB数据库与华为DPA数据保护一体机兼容性互认证。 图1：华为DPA数据保护一体机兼容性互认证 亚信安慧AntDB数据库作为领先的数据库解决方案提供商，专注于数据库产品的研发与创新，以其卓越的性能和稳定性，服务于超数亿用户，连续十年无故障运行。亚信安慧AntDB数据库的云原生分布式架
阅读更多...

Copyright China编程 23002807@qq.com

沪ICP备2023033072号-2