『TBtools-plugin』OneStepWGCNA插件

这个插件已经下架了,上架了更好用的WGCNA shiny:https://www.jianshu.com/p/4691c535a522

常见问题反馈和总结

插件释放了有差不多一个月了吧,总结了下大部分同学反馈回的问题大致分为以下几种情况:
以后类似问题不做回答。
先放一下示例数据 感谢王瑞老铁提供。如果用示例数据跑出结果了,说明插件安装完成。具体设置如下图

20210305131252.jpg

提示成功后输出结果文件夹下有以下文件说明没问题了。


20210305130349.jpg
  • 没有装 Rserver插件! 会直接报错,Rserver相对于在TBtools中装了个R语言,而基于R语言写的插件运行也就是在TBtools的R语言环境中进行Rscript xxx.R如果你连R都没装,怎么可能运行起来? Rserver插件的获得方法目前只有CJ在TBtools用户群中释放。
  • 包安装的问题: 如果在test log中发现报错信息中包括下面信息的就是包没有装成功,这个大概率是网络问题,换个网络再运行几次估计就可以了。比如你可以连接你的手机热点。反正我连河南移动基本没问题。
Error: xxx...package....
there is no package called 'pillar'
  • count data和trait data有NA值或者空白: 如果你编辑traitdata的时候多打了两行空白行,或者在某行后面多安了个tab,都会造成外表上看见你的trait data和count没有问题,实际上读近R之后这些空白还是会被读入data.frame中的,肯定会报错,不会撸代码的话就人肉排查把。把那些乱七八糟的都删掉...

  • 你想保留的基因数量超过了过滤后的基因数量:
    如果你仔细看了下面博客每个参数详细的介绍,这是很好猜到的,比如你的表达量矩阵中共有5w个基因(这里忍不住提醒一下,建议大家的表达量矩阵用所有的基因,不要挑差异基因或者自己选基因做,具体原因自行查文献,看WGCNA操作手册,FAQ,google ) 我写的流程中会进行第一轮筛选,就是下面说的在几成样本中fpkm或者readcount大于n,如果你第一步设置的阈值太高,比如九成样本中FPKM > 100,这样过滤完以后可能就剩下3000个,然后在第二步过滤的时候你选择保留8000个,就会报错。这个时候适当的减少阈值,比如使用readcount的时候用5成样本中readcount > 10或者5成样本中FPKM > 1这些,当然这里是可变的,可以自己调整。

  • 最后一种我想骂人,拜托你们也尊重下作者,好好看看下面的文章!
    有几个老铁给我邮件说跑不通,报错整的我也没看明白,结果浪费了我好久最后发现他好几个参数都没有填!!! 没有选参数怎么能跑出来!所以如果没有仔细阅读下面文章,特别是参数介绍(脑图那块)出现报错的,麻烦您再回去看一遍!!!


国际惯例先上图


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还是fpkmreadcount是不可以直接用来做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 cutoffSamplePerc: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的机会...

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 212,332评论 6 493
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,508评论 3 385
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 157,812评论 0 348
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,607评论 1 284
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 65,728评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,919评论 1 290
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,071评论 3 410
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,802评论 0 268
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,256评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,576评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,712评论 1 341
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,389评论 4 332
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,032评论 3 316
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,798评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,026评论 1 266
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,473评论 2 360
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,606评论 2 350

推荐阅读更多精彩内容