『TBtools-plugin』OneStepWGCNA插件
国际惯例先上图
OneStepWGCNA
再放链接：
Gitee repo - ShawnMagic / ShawnTBPluginR
写在前面的一些废话
刚从北京上完课回所那会，何老师建议我学下WGCNA，可能对我可课题有用，当时R还没怎么入门（还逃过一次赖江山老师的课，哈哈）于是就准备开始鼓捣WGCNA，说实话，那时候学习软件全靠简书，知乎，CSDN几位大佬的教程，估计和大多数朋友一样，入门的都是陈同老师的这篇WGCNA分析，简单全面的最新教程,那时候也没有看官方文档的习惯，说实话，连里面的术语啥的都不懂，什么加权共表达基因网络，什么SFT，什么KME... 当时就想着follow大佬的代码跑通是正事，后续Jimmy大佬的一个直播他自己学习hope大佬的教程视频 我是如何学习WGCNA分析，对我启发很大，看文献，看manual，试代码....
我自己比较笨，大概花了半年才逐渐的理解了所有的代码... 但由于数学硬伤，统计，图论那一块真的搞不明白.. WGCNA的原理到现在感觉还是在云雾中，直接看文献基本等于在看天书... 可能需要毕业后静下心来，好好攻克一下。
之后的事就简单了，前后做了好多次，随着潜移默化的发展，R写代码的水平回看起来比之前有了不小的提升，所以后来把分析切分成了几个大模块，也就形成了今天插件的雏形：
- 数据清洗部分
- SFT，power计算，样本聚类
- 一步法网络构建
- module-trait关系
后续的分析基本上就是拿出感兴趣的模块，找hubgene，导出关系，富集分析等等。那些个性化的操作不想一股脑把所有的模块都搞出来，主要是占地方，一下出一串文件大部分也没有什么卵用。所以后面会写新的插件实现这些功能，根据前面找出的感兴趣的模块，后续再进行分析。
插件说明
首先说明一个问题，流程化的东西，多少都会有弊端，所以这个插件只适合粗略的进行WGCNA分析，大致看一下模块，trait-module的关系。因为WGCNA分析确实不是那种一键就能搞出非常牛逼的分析结果... 因为有太多的参数需要调整了。如果想要真正的做出一个漂亮的WGCNA结果。
建议必须阅读以下文章：
Tutorials for the WGCNA package
WGCNA-FAQ 可以解决你80%的问题
WGCNA: an R package for weighted correlation network analysis
当然还有简书，知乎各位大佬的文章，这个就靠你自己去搜索了。没有好坏，各有特点！
插件代码主要follow了一下两位大佬的代码,如果想快速了解建议阅读：
WGCNA分析，简单全面的最新教程
STEP6:WGCNA相关性分析
使用方法
首先按照CJ大佬的教程安装Rserver插件这个神级插件是所有基于R的TBtools插件的爸爸😄，相当于在TBtools里面装了个R，
然后在我的Gitee仓库下载OneStepWGCNA插件这个后续有更新的话我会及时推送，欢迎各位反馈bug.... 课题组的老铁给我说了俩..已经改正...
和TBtools使用的规则一致，拖进去，选参数，点点点...
首次运行时间较长，应为要装一堆包，包装不上大概率是因为网络！实在不行就连热点！
参数详解
OneStepWGCNA
OneStepWGCNA
主要有7个参数：
-
readcount or fpkm这个文件是表达量矩阵，不能有多于的信息，列名（column name）是sample名称，行名是geneID，中间的数据是FPKM值或者readcount 如下表
| GeneID | Sample1 | Sample2 | Sample3 | Sample4 | Sample... |
|---|---|---|---|---|---|
| Gene1 | 1 | 2 | 1 | 2 | ... |
| Gene2 | 0 | 1 | 0 | 2 | ... |
| Gene3 | 2 | 3 | 2 | 1 | ... |
| Gene... | ... | ... | ... | ... | ... |
-
Trait data这里分两种，一种是数量性状（比如百粒重，株高，生育期等这种可以计数的）还有一种是质量性状（种皮颜色，有无表皮毛，有无腺体等），或者分类，第一种数量性状表列名为表型，行名为sample name，第二种只有两列第一列为sample name，第二列为分类特征。如果既有数量性状，又存在质量性状，建议质量性状的有无用0和1代替（脚本自动处理也是根据分类生成0-1矩阵）
数量性状表
| GeneID | Trait1 | Trait2 | Trait3 | Trait4 | Trait... |
|---|---|---|---|---|---|
| Sample1 | 1 | 2 | 1 | 2 | ... |
| Sample2 | 0 | 1 | 0 | 2 | ... |
| Sample3 | 2 | 3 | 2 | 1 | ... |
| Sample... | ... | ... | ... | ... | ... |
分类表
| Sample | Type |
|---|---|
| Sample1 | case |
| Sample2 | case |
| Sample3 | case |
| Sample4 | control |
| Sample5 | control |
| Sample6 | control |
| Sample... | … |
-
Expression data这里意思是你输入的表达量矩阵到底是readcount还是fpkm，readcount是不可以直接用来做WGCNA的，需要标准化，根据官方和Jimmy大佬学徒的帖子我测试后推荐了两个选项,一个是基于DESeq2的varianceStabilizingTransformation, 另一个是基于edgeR计算的log10(CPM+1)有老铁反馈lgCPM的结果不太好，说看看log2(CPM+1)怎么样,看大家的反馈结果... 如果用FPKM/RPKM/TPM这类已经标准化后的值，这里给了两个选项，一个是做log另一个是原始的值，这里众说纷纭，有的说log后把很多差异都抹掉了，还有的说不取log会放大差异...我也是头疼.. 到底哪种好，我自己感觉你都试试... 能搞到合适的power，构建一个真正的scale-free网络他就是好方法。😂 -
RcCuttoff 和 SamplePerc在WGCNA-FAQ中作者提到，最好去除一下noise，再去做，（Jimmy大佬也说这不太重要，重要的还是你选多少个基因进行WGCNA）官网建议的是保留九成以上样品中readcount>10的基因，当然后面还有强调就是根据测序深度，自己的生物学问题而改变。那么这里就给出了两个参数RcCutoff： readcount cutoff和SamplePerc：sample percentage,之前我都是基于Readcount做的，所以当时取变量名的时候顺手就写成这个了，估计不会改了，fpkm/RKPM/TPM也是适用的，比如你想5成样品中fpkm大于0.2那么RcCutoff那里就填0.2，SamplePerc那里就填0.5。 -
RemainGeneNum:这个参数就是你需要保留多少个基因进行WGCNA，WGCNA并不是越多越好，当然你基因足够大，你用的PC也很难全部跑完，对内存的要求比较严格，跑得多了你会发现当基因数量到达一定阈值之后，再增加只会使grey模块变得更大，因为在没有样品异质性的RNAseq中，大部分基因的表达都是保守的，意思就是基本在所有样品中他们的表达量是一致的。增加再多的也没有太大用处。这里采取的筛选方法还是鲁棒性较强的MAD，先算MAD，再rank，选取MAD大于0.01的前RemainGeneNum个基因作为WGCNA的input。这里有个联动各位注意！！ 前面在去噪音的时候会去掉一部分基因，如果你选择条件太苛刻，可能会去掉一大半，比如昨天代一个师弟做，他选择9成的FPKM大于10，7w个基因去的只剩下几千个...，后面这里又设置了1w个，所以一炮就报错。所以一定要注意
结果文件
-
01.test.Sample_clustering 样品聚类结果
20210117114017.jpg
-
这里如果出现明显的离群样本，记住样本的名字，在表达矩阵和trait矩阵中删掉对应样本再跑一次
-
01.test.SFTPlot.pdf 软阈值筛选
20210117114055.jpg
-
绿线为0.8 红线为 0.9，脚本内筛选值为0.8
-
02.test.Eigeng_adja_heatmap 邻接基因矩阵热图
20210117114126.jpg
-
反映了模块相关性（模块小于等于3没法做）
-
02.test.Module 模块划分
20210117114156.jpg
- 03.test.Module_trait 模块-矩阵热图
20210117114233.jpg
- 还有几个Rdata文件，是保存了几个变量，在后续的WGCNA-module分析插件中有用，所以现在先保存了一下，倒是后直接load
- test.log保存的日志，脚本运行的日志，如果保存可以查看原因。
- 两个表格，一个是module-gene对应关系，另一个是所有基因的连通性
最后再唠叨两句
做WGCNA最主要的是一定要清楚WGCNA是干什么用的，第二是你的科学问题是什么，常见的困扰大家的问题有几个我总结了下：
- 到底是readcount还是fpkm，这里见仁见智，还是我说的，你搞出robustness高的scale-free网络来是终极目标... 没人去纠结到底input哪个更科学，Stack Overflow上有大佬直接用z-score的fpkm也做出来了....所以input方法千变万化，只要有道理就可以，不可死板教条，要灵活
- 搞不到合适的power.. 这时去查看下你的input，有很多时候大家不清楚WGCNA的作用，一股脑把所有的表达量都扔进去，没有合适的power大概率是由于样品异质性（heterogeneous）造成的，或者换句话说你把不同组织的样本混合一起跑，其实在同一物种内组织间的差异是要远远大于品种间的差异的，你这么做如果生物学问题是想看不同组织间特异表达的基因模块没有问题，如果你想看不同材料间的，建议洗洗睡吧（如果样本够建议分开组织各自做）... 所以一定要注意自己关注的生物学问题，减少变量。
- module-trait相关性低，很好奇文章里姨妈红是怎么来的。我的经验... 如果直接拿数量性状关联，大概率会相关性不高。我的做法是把数量性状分级，整体做个数据分布图，设置一下bin的大小，按bin分级，或者boxplot按照分位数分级。把数量性状转换成等级... 效果确实会好一点，这个也好理解，如果你做的品种材料间的性状，不是野生型突变体之间，或者近等基因系之间的，差异不会太显著，波动范围很小，基因的表达量的波动也很难和你性状数值的波动相吻合，但是你分级之后应该会变好。只是我个人的经验不一定对。
- hubgene怎么搞，目前方法还是众说纷纭，这部分会在下一个插件中推荐一个更好的方法，但是！ 一定要注意前提，就是你在前面的分析中得到了一个合适的power，构建了一个scale-free网络，有些分析搞hubgene是没有意义的，比如强烈样品异质性的WGCNA分析，划分出module就到头了，应为你的网络很难成为scale-free网络，那么hubgene也无从谈起了，你会发现模块内基因的连通性大部分都在0.9以上... 原因就是你再做WGCNA之前已经把模块划分好了，基本上是一个组织对应一个模块...
综上TBtools牛批，又给我一个装x的机会...
