前言
Smin是在蛋白质功能预测中比较流行的一个指标,具体由来我也不甚清楚,只是在最近复现的几篇论文中反复出现了,所以记录一下。
计算方法
(图来自于PSPGO论文)
其中𝜏表示阈值,t表示GO标签,Pa(t)表示该GO标签依赖父亲。
Prob(t|Pa(t))使用频率估算,是描述整个Protein-GO数据集的统计值。
Pi(𝜏)表示在阈值𝜏下的第i个蛋白质的预测GO集合,Ti表示第i个蛋白质实际GO标签集合。
这个计算流程比较复杂,得先用频率估计概率求出ic(相信大部分GO预测论文里面都有对应代码)
然后枚举Threshold值,对每个蛋白质计算 当前预测GO标签传递闭包 与 真实标签的传递闭包 的差集(Ti-Pi的差集和Pi-Ti的差集),统计这两部分的ic值各自的平均数,然后求这两个平均数的评分平均数作为S值。
最后找到最小的S值就是Smin了。
一些疑问的解答
1、为什么求传递闭包?
(什么是传递闭包?请看《离散数学》图论部分)
首先GO标签依赖关系是一个DAG(这里只考虑is_a和part_of关系,加入regulates之后就可能存在环了!)
然而一般数据库给出go.gaf文件都只是给出最具体的标签部分,更上层的GO标签未必有记录,但是有子标签必然是有父标签的,所以我们需要读取go.obo文件建立DAG图,然后求每个子标签的传递闭包的并,这才能得到一个蛋白质完整的GO标注数据。
2、ic值的计算
首先ic值也必须建立在所有蛋白质的GO标注都求过传递闭包的基础上。
在一个点有两个父亲的情况应该怎么办,例如:
在一些论文中,我们选取的是Prob(t|Pa(t))更大的作为这里标签的ic值(也就是选取父亲出现次数更少的进行计算)。另一些论文中可能会有不同的处理。
3、计算Smin的时候能否先对预测结果做一遍传递闭包?
复现论文是发现的细节,有些论文在计算Smin的时候巨慢无比,原来是每次都会把预测结果传递一遍再来求ru和mi值。
这样按照理论,得到的Smin值会更优一些,因为可以避免掉一些预测误差导致的中间某个GO标签缺失的情况。
实战来说,确实可以这么做,毕竟预测结果肯定是可以先传递一遍,再给别人拿去用的。
但有些老实的论文并没有做传递,得到的Smin值可能就会略有偏高。
4、关于Smin计算的优化
在复现PO2Vec论文时发现的问题,算一个thresh得花我十多分钟,实在受不了了,就给他优化了一把。
优化思路:提前保存每个GO标签的传递闭包结果,更改计算顺序,利用two-pointers计算每个蛋白质对各个阈值下的ru与mi的贡献。
在优化前,在swissprot数据集下,计算一次mf的Smin需要约10个小时,在优化后,只需要不到一个半小时。
主要优化的evaluate_cafa3.py中的evaluate_model_predicion函数。
def evaluate_model_prediction(labels, terms, model_preds, go_rels, ont):ru_list = []mi_list = list()# go setgo_set = go_rels.get_namespace_terms(NAMESPACES[ont])go_set.remove(FUNC_DICT[ont])# labelslabels = list(map(lambda x: set(filter(lambda y: y in go_set and y in terms, x)), labels))ancs = {}for go_id in tqdm(go_rels.ont,desc="calc ancs..."):ancs[go_id] = go_rels.get_ancestors(go_id)lim_set = go_set.intersection(terms)ru_list = [0]*101mi_list = [0]*101total = 0for x in labels:if len(x) > 0:total+=1assert len(labels) == len(model_preds)for label,pred_score in tqdm(zip(labels,model_preds)):tmp_lis = list(zip(pred_score,terms))tmp_lis.sort(reverse=True)j=0new_annots = set()for t in reversed(range(101)):threshold = t/100.0while j<len(tmp_lis) and tmp_lis[j][0]>threshold:new_annots |= ancs[tmp_lis[j][1]]j+=1new_annots = new_annots.intersection(lim_set)tp = new_annots.intersection(set(label))fp = new_annots - tpfn = label - tpfor go_id in fp:mi_list[t] += go_rels.get_ic(go_id)for go_id in fn:ru_list[t] += go_rels.get_ic(go_id)ru = np.array(ru_list)/totalmi = np.array(mi_list)/totalsmin = np.min(np.sqrt(ru * ru + mi * mi))return smin