GO Level 到底是什么?
最近忙着鼓捣硕士毕业论文。。。上周一重新确定了论文题目《XXX数据库的构建与应用》...于是用了大半周的时间搭了数据库和网站,然后开始整毕业论文。于是今天才上来看公众号,查看后台,连续看到两则非常有趣的信息,贴出如下:
长这么大,我最担心的,就是被别人看穿。今次,这两个事情都很有趣,两个朋友先后看穿了我:
1. 我确实是认真的看过了健明鼓捣的论坛。此外,借此机会,我还要再声明一次。我觉得应该算是我们小群里面,不算唯一支持,也算是最支持这位朋友提到的论坛的,这个在之前健明托了我们几个人微信小聊的时候我就明确表示过啦。后来在omicshare论坛线下交流,我甚至还推荐了这个论坛。但是,由于我们喜欢的群体不同,所以我也就没进去鼓捣了。
2. 我确实不愿意分享。主要原因,我发现随着年龄的增长,我的工作效率越来越低,该做的想做的依然忙的不可开交,自然不可能有时间分享。其次,很久很久以前,我写过一些博客,也写过一些小工具。我还记得早前写得一些博文,比如GOstats做GO富集和KEGG富集,pheatmap的使用等软件使用文档,或者是转录组数据分析流程等文档,在我博客关闭很久很久之后,依然有人找我交流,我基本都还是正常回复的。我不愿意分享,因为确实不知道分享是什么。总之这个分享,不会是写博文,也不会是专心鼓捣好一个QQ群,也不会是经常混迹在论坛与人交流,更不会是写TBtools。
===========以下贴上很久以前的博文=================
What is GO term level
经常会在一些文献上看到so-called GO term level3/level2 的柱形图绘图
以前写过一个blast3go的小玩具,可以实现GO idmapping 从NCBI的gene id映射为GO id
进行so-called level2/level3 绘图,以及常用的GOstat topGO等富集分析,出图
还是有几个人比较喜欢这个玩具的,虽然,因为有大神已经写出了超级流弊的工具,我决定放弃往哪个方向继续玩了
感兴趣的同学可加入我们的QQ群来交流:bioinformatics*中国(276151571)
===================
但是回到正题,blast3go,其中实现level2/level3 的绘图部分,使用的是非常机械的方法
具体做法如下
1. 使用bioconductor的GO.db这个强大有及时跟新的GO包获取MF BP CC三个term的子term
2. 将注释数据导入,packageAllGOFrame, 导出
3. 使用1.获得的level2/level3列表,对2.导出的数据进行计数统计,使用ggplot2绘图
以上做法,可以接受,但是毕竟还是不够流弊,anyway,起码比一直没更新的WEGO的全面
由于项目需要,昨晚重新查看了’到底GO level’是个什么东西
GO 官网是最权威的
GO terms do not occupy strict fixed levels in the hierarchy. Because GO isstructured as a graph, terms would appear at different 'levels’ if differentpaths were followed through the graph. This is especially true if one mixes thedifferent relations used to connect terms. Thus it is more proper to ask: “whatis the maximum depth of such and such a term” (or minimum, average, etc.).
GO官网的说法,根本不存在所谓的GO level,,,
但是为什么blast2go上面有GO level分布信息?怎么算的
其实了解GO ontology(我习惯称为层级结构)的同学应该很清楚,所谓level ,既然官网这么说了,那么要么就是最长路径要么就是最短路径或者就是平均路径
我们再观察Amigo对GO term ontology的摆放,第一层,第二层,第三层,,,
# 从刚提取出来的max depth信息
# GO:2001275positive regulation of glucose import in response to insulin stimulusbiological_process 8 # 第8层
非常直观地发现,应该就是最长路径, 图片上,明显按照最长路径(max depth)先排布,随后使用另外的箭头标出捷径
于是我们就从本地化的mysql数据库中,直接导出所有GO term 的最长路径(我们可以认为这就是官网语句提到的,maximum depth)
$ mysql -ub2gUserName -pb2gPassword
mysql>/
SELECT DISTINCT
term.acc,
term.name,
term.term_type,
max(p.distance)
into outfile '/tmp/b2gMaxDis.txt'
FROM
term
INNER JOIN graph_path AS p ON (p.term2_id=term.id)
INNER JOIN term AS root ON (p.term1_id=root.id)
WHERE
root.is_root=1
AND term.is_obsolete=0
GROUP BY
term.acc,
term.name,
term.term_type;
mysql>exit;
$ mv /tmp/b2gMaxDis.txt.
$ perl -ane '@F=split /\t/,$_;print $F[3]' b2gMaxDis.txt|sort|uniq -c
————–
1 0
3 1
1808 10
935 11
316 12
109 13
8 14
59 2
341 3
2046 4
4788 5
9541 6
7660 7
8796 8
4295 9
————-
就是只有14个level…与blast2go的结果一样,就是到第14个level
我们反观这个做法,对照下我们blast3go的做法 …
perl -ne '@F=split /\t/,$_;print $F[2],"\n" if $F[3]==2' b2gMaxDis.txt
blast3go还多出了9个,难道是我mysql数据库就没更新过来?
anyway,我还是坚持blast3go的做法,毕竟,,,反正我觉得bp mf cc的直接子term就是最值得关注的,管你是不是全都都是所谓的level2/level3
这个选择关系到统计出来的信息是否有效,如果我们就按照max depth选择了level 3, 这个可想而知,有些term,虽然是level 2/3
同个level 对应的信息量,本来不等,虽然bp mf cc对应的这9个子term被包含了其他信息,但其还是足够具体…. Whatever, who care?