以下文章来源于生信菜鸟团 ,作者Quasimodo
生信菜鸟团荣誉归来,让所有想分析生物信息学数据的小伙伴找到归属,你值得拥有!
其实转录组走到现在我总觉得少了点什么东西,后来才想起来是cytospace寻找hub基因
这个其实很简单,网上的教程大把大把的,而且不需要写代码
参考:
RNA-seq入门实战(十):PPI蛋白互作网络构建(下)——Cytoscape软件的使用
Cytoscape教程|如何挖掘蛋白互作网络中的核心基因?(https://zhuanlan.zhihu.com/p/595214085)
寻找核心基因+子网络 (https://zhuanlan.zhihu.com/p/137653031)
但是随着我学习到的一些其他知识,让我加深了这方面的理解,所以还是决定写一下
所以这并不是一篇Cytoscape教程,而是一片探索性的文章,我也没有把推文标题写成寻找PPI中的hub基因,有需要的同学可以参考上面给出的其他内容
我这里一开始是用的是 初探mRNA、lncRNA联合分析之下游 的差异表达lncRNA共表达mRNA对应的基因
打算是进一步narrow down :一是看看PPI网络和共表达网络有什么不同,二是找出来的PPI hub gene又是哪些
差异表达lncRNA共表达mRNA:
SWT1
AXIN2
HARS2
MMS19
HYLS1
PTRH1
MXRA7
PUF60
ASAH1
CLEC16A
NFIX
ZNF638
FIGNL1
BID
CNOT2
PTPN12
CD44
NCAM1
CANT1
DUSP14
REPS1
MFAP3L
EPB41
TCF4
IFNGR1
RNASEH2B
之所以说是进一步,是因为我们一般拿到共表达网络就可以了,不会再对其进一步做PPI找hub
在这篇推文中我们已经进行了
并做了比较
直接string检索可以发现有5个基因存在关联,特别是CD44处于中心位置,这和前面的lncRNA共表达mRNA的结果是一致的,可以发现CD44的连线很多
此外,还有EPB41 TCF4 AXIN2 NCAM1 ,其中EPB41 TCF4属于前面WGCNA、lncRNA共表达mRNA的overlap。而实际上前面WGCNA、lncRNA共表达mRNA的overlap的基因位置还比较特殊,属于右下角小局部网络中,这也进一步提示我们关注这个共表达网络中的小局部网络(不光WGCNA和lncRNA共表达mRNA的overlap集中在这里,而且所有共表达的mRNA的PPI hub中5个有两个位于这里,这同时反映了表达谱特征和真实生物网络)。
WGCNA、lncRNA共表达mRNA的overlap(黄框为overlap,图为共表达网络):
有意思的是,AXIN2 NCAM1在string中的PPI关联是通过文本挖掘找到的,而并没有找到共表达信息,而在前面的lncRNA共表达mRNA网络中,这些基因(CD44 AXIN2 NCAM1)是通过了作者设置的与差异表达lncRNAs共表达(转录水平上)的阈值
The absolute values of correlation coefficients >0.95 and p < 0.001 were specified as screening criteria.
若导到cytospace中
#node1 node2 node1_string_id node2_string_id neighborhood_on_chromosome gene_fusion phylogenetic_cooccurrence homology coexpression experimentally_determined_interaction database_annotated automated_textmining combined_score
AXIN2 CD44 9606.ENSP00000302625 9606.ENSP00000398632 0 0 0 0 0 0 0 0.546 0.546
CD44 EPB41 9606.ENSP00000398632 9606.ENSP00000345259 0 0 0 0 0.042 0 0 0.630 0.630
CD44 TCF4 9606.ENSP00000398632 9606.ENSP00000381382 0 0 0 0 0.104 0 0 0.527 0.558
CD44 NCAM1 9606.ENSP00000398632 9606.ENSP00000480132 0 0 0 0 0 0 0 0.712 0.712
太少了
拿参考文献提供的差异表达的mRNA做,顺便看看这样PPI cytospace 寻找hub基因的结果和作者基因水平的有哪里不一样
https://ftp.ncbi.nlm.nih.gov/geo/series/GSE209nnn/GSE209825/suppl/GSE209825_mRNA_ADS_vs_Normal.Differential_analysis_results.xls.gz
The absolute logarithmic fold change (LFC) ≥ 1 (|log2 fold change|≥ 1) and p-value <0.05 were adopted as the screening criteria
55 awk '$(NF-1) < 0.05 && $(NF-2) < -1 {print $3}' GSE209825_mRNA_ADS_vs_Normal.Differential_analysis_results.xls | sed '1d' > DEG_list.txt
56 awk '$(NF-1) < 0.05 && $(NF-2) > 1 {print $3}' GSE209825_mRNA_ADS_vs_Normal.Differential_analysis_results.xls | sed '1d' >> DEG_list.txt
$ wc -l DEG_list_uniq.txt
2492 DEG_list_uniq.txt
超过2000,直接导入cytospace
使用MCODE寻找子网络(模块):
原文献:
In order to construct a visual network map, interaction relationships for the list of differential genes were extracted from the STRING protein interaction database (https://www.string-db.org/). The network data files were imported directly into Cytoscape software for visual editing. The size of a node in a PPI map is proportional to the degree of the node. Among them, PPI nodes with relatively high connectivity include UBA52, AKT1, SUPT20H, RPL19, EGF, and MYC
relatively high connectivity 指的是什么算法呢?
看起来有点像最简单的degree
反正作者应该是直接用对整个网络找的hub,没在这之前划分子模块,所以我们后面还是对整个网络进行hub的寻找
cytoHubba: identifying hub objects and sub-networks from complex interactome (https://bmcsystbiol.biomedcentral.com/articles/10.1186/1752-0509-8-S4-S11)
cytoHubba提供了12种算法,虽然我在其一开始发表的文章中只找到11种,这里还是归纳一下:
在网络分析中,径向度是一种用于衡量节点在网络中的中心性的指标。它关注的是节点与其可达邻居之间的距离关系。具体来说,径向度将高中心性赋予那些在其可达邻域中离其他节点的距离较短的节点。
为了计算节点的径向度,首先需要确定节点的可达邻域。可达邻域包括节点本身以及所有与该节点距离不超过其直径的节点。直径是指可达邻域中节点之间的最大距离。
然后,对于每个节点,计算它与可达邻域内所有其他节点之间的平均距离。距离越短,节点的径向度越高。因此,径向度高的节点可以被认为是网络中的中心节点,它们在网络的结构和连接中起着重要的作用。
通过计算节点的径向度,我们可以获取关于网络中节点的重要性和中心性的信息,这对于理解网络结构和节点的功能扮演了重要的角色。
可以将瓶颈定义为具有高介数中心性的蛋白质(即,具有许多“最短路径”的网络节点,类似于高速公路地图上的主要桥梁和隧道)。
参考:The Importance of Bottlenecks in Protein Networks - NCBI (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1853125/)
该方法基于网络在信息流传递过程中的瓶颈原理,寻找在信息流中具有关键作用的基因。Bottleneck方法关注的是基因在网络中连接不同模块之间的重要性,例如,一个基因连接了多个模块,其在信息传递中起到了关键的桥梁作用。因此,Bottleneck方法寻找的是在网络中起到关键枢纽作用的基因。
上面带了引用没给参考的部分都是AI帮我解释但我不确定的,需要甄别着看,也希望懂的老师能在评论区说说,或者笔记投稿,我在检索上面这11个方法的时候其实无论中外互联网感觉相关资料都很少
挑选一些方法寻找top10的hub nodes:Radiality Bottleneck stree betweenness EPC degree MCC ACTB ACTB ACTB ACTB ACTB ACTB RPL15 AKT1 AKT1 AKT1 AKT1 AKT1 AKT1 RPL17 CDC42 CDK2 CTNNB1 CDC42 GAPDH CTNNB1 RPL19 CTNNB1 CTNNB1 GAPDH CTNNB1 HSP90AB1 GAPDH RPL6 GAPDH FYN HSP90AB1 GAPDH MYC HSP90AB1 RPLP0 HSP90AB1 GAPDH MYC MYC RACK1 MYC RPS16 MYC MYC RPS27A RPS27A RPS27A RPS27A RPS2 RPS27A RPS27A TRIM28 TRIM28 UBA52 UBA52 RPS20 UBA52 TRIM28 UBA52 UBA52 UBB UBB RPS3A UBC UBC UBC UBC UBC UBC RPS9
其实可以发现全局算法的结果较为一致,最简单的局部算法degree也差不多,但是MCC我们发现找到的都是核糖体大小亚基蛋白
随着调整不同算法出图,我发现作者的PPI hub中始终有一些没有出现在我的top10中
检查:
> PPI_hub<-c("MYC","EGF","AKT1","SUPT20H","UBA52","RPL19")
> PPI_hub %in% deg_list_uniq
[1] TRUE FALSE TRUE FALSE TRUE TRUE
可以发现EGF SUPT20H在转录水平的表达差异并不显著,而在作者的DNA水平定量上显著,所以我们构建的PPI网络中没有这两个节点。其实,这也在侧面说明了在真实生物世界中网络调控处于中心位置的蛋白并不一定差异表达非常显著,可能仅仅是略微的差别就有可能引起一系列变化从而影响表型。
在检索学习如何鉴定网络中心节点的时候,我检索到了一个整理的很好的网页,可以对上面11种方法补充着看,这里也贴出来
网络中心性 (https://wiki.swarma.org/index.php/%E7%BD%91%E7%BB%9C%E4%B8%AD%E5%BF%83%E6%80%A7#.E6.B8.97.E6.BB.A4.E4.B8.AD.E5.BF.83.E6.80.A7_Percolation_centrality.EF.BC.88PC.EF.BC.89)
在图论和网络分析中,中心性 centrality 指标用于识别图中最重要的顶点。其应用包括在社交网络中识别出最有影响力的个人,在因特网或城市网络中识别出最为关键的基础设施节点,以及识别疾病的超级传播者。中心性的概念最初是在社交网络分析中发展起来的,许多用于衡量中心性的术语都反映出了它们的社会学起源。中心性不应与节点影响度相混淆,后者意在量化网络中每个节点的影响。
其中包括了一些其它中心性,如调和中心性、特征向量中心性、网页排名中心性等
其中,网页排名中心性 PageRank centrality
PageRank算法是搜索引擎Google所采用的评估网页重要性的算法,GeneRank算法在PageRank的基础上通过生物网络的拓扑结构和节点的初始值对节点分数进行迭代,可以成功地对网络中的基因进行优先排序
original GeneRank:
Modified GeneRank:
我们这里不做公式的深入解析,只是借此引入节点的初始值
虽然我们在前面发现
EGF SUPT20H在转录水平的表达差异并不显著,而在作者的DNA水平定量上显著,所以我们构建的PPI网络中没有这两个节点。其实,这也在侧面说明了在真实生物世界中网络调控处于中心位置的蛋白并不一定差异表达非常显著,可能仅仅是略微的差别就有可能引起一系列变化从而影响表型。
但是初始分数(log2FC的绝对值或者-lgPvalue)其实还是能在一定程度上说明节点的在真实的生物世界网络中的重要性
私以为其实可以两个都做做,先不带初始值看看hub节点是哪些,这时的hub节点都是因为在网络中占据了重要地位。然后带上初始值,就可以看看哪些节点是因为带上了初始值而变成hub
思路打开
另外这些无初始值的算法并不一定需要用cytospace不可,R包 centiserve已经整合了
所以不一定非要是PPI,只要是网络结构就可以做
思路打开
现在来看看有初始值的PPI如何找到hub节点
string最开始是没有这个功能的,后面给加上了
但是这个功能和我一开始想的并不一样
我以为的是可以加上初始值来确定hub基因了
其实只是一个类似GSEA的功能,关于基因集分析可以参考我前面的这篇系统介绍的推文基因集分析的前世今生(附进行通路富集分析的9个tips)
测试数据结果:
看来string和cytospace还是分割的,功能联系没那么紧密啊
我手上有21个基因及其初始值的列表输进去看看
“具有值/等级的蛋白质”功能 基于大基因列表(覆盖整个基因组或至少其大部分)中的值的分布进行统计富集检测。
考虑到我们进行的测试数量,几乎不可能在如此少量的值上产生具有统计意义的信号。然而,您仍然可以继续并将您的数据集提交给基于常规基因集的分析,我们将您的蛋白质集与整个蛋白质组进行比较
外面的环就代表了初始分数
相当于没用初始分数,只是让你看看在对应的网络中哪些节点处于中心位置并且初始分高(倒也是一种思路)
但我还想看看初始分数最后对hub基因的影响
想要找hub基因还是需要导出来去cytospace,但是这个过程就丢失掉了nodes初始分数,只是edges发挥作用
我们刚刚介绍的GeneRank也是需要输入初始值的
我们可以自己手搓一下上面那个GeneRank算法,并不难
鉴于R的运行速度,还是推荐用相对来说快一点的python
如果熟练可以自己用基本数据结结构比如列表、字典完成整个任务,此外我发现NetworkX这个第三方库用起来还可以,如果用python画PPI可以考虑一下这个
https://networkx.org/
总之,在完成算法迭代前,我的初始分数top10为
CCR1, IFNAR2, MAPT, OAS1, HLA-C, ABO, NSF, PDE4A, CAT, FUBP1
完成算法迭代后top10为
max_dict = {}
for index, value in enumerate(ls_b):
max_dict[value] = index
sorted_dict = sorted(max_dict.items(), reverse=True)[:10]
top_10_indexes = [index for _, index in sorted_dict]
print(top_10_indexes)
[13627, 3767, 1609, 2927, 6738, 482, 452, 1953, 11160, 3862]
selected_elements = [nodes[index] for index in top_10_indexes]
print(selected_elements)
['ABO', 'CCR1', 'IFNAR2', 'MAPT', 'OAS1', 'HLA-C', 'JAK1', 'NSF', 'CYP4B1', 'PDE4A']
可以发现,网络的拓扑结构特征被融入到了基因中,特别是JAK1和CYP4B1这些一开始并不在top10的节点
这样的top节点作为hub基因会不会更有意义呢?
那是那些句话,这里就刚好放到最后作为小结:
EGF SUPT20H在转录水平的表达差异并不显著,而在作者的DNA水平定量上显著,所以我们构建的PPI网络中没有这两个节点。其实,这也在侧面说明了在真实生物世界中网络调控处于中心位置的蛋白并不一定差异表达非常显著,可能仅仅是略微的差别就有可能引起一系列变化从而影响表型。
但是初始分数(log2FC的绝对值或者-lgPvalue)其实还是能在一定程度上说明节点的在真实的生物世界网络中的重要性
私以为其实可以两个都做做,先不带初始值看看hub节点是哪些,这时的hub节点都是因为在网络中占据了重要地位。然后带上初始值,就可以看看哪些节点是因为带上了初始值而变成hub
无论是WGCNA、PPI+cytospace的各种方法、调控序列共表达网络还是其它考虑了初始权值的方法,其实都可以根据自己的实际情况来,思路打开没有坏处,narrow down将会是高通量数据的主旋律,就算是学到后面用机器学习(回归模型、随机森林、支持向量机等等)也是同样的思路
微信扫一扫
关注该公众号