依赖标签分类任务Smin值计算（蛋白质功能预测，GO标签）

前言

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 set
    go_set = go_rels.get_namespace_terms(NAMESPACES[ont])
    go_set.remove(FUNC_DICT[ont])
    # labels
    labels = 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]*101
    mi_list = [0]*101
    total = 0
    for x in labels:
        if len(x) > 0:
            total+=1
    assert 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=0
        new_annots = set()
        for t in reversed(range(101)):
            threshold = t/100.0
            while j<len(tmp_lis) and tmp_lis[j][0]>threshold:
                new_annots |= ancs[tmp_lis[j][1]]
                j+=1
            new_annots = new_annots.intersection(lim_set)
            tp = new_annots.intersection(set(label))
            fp = new_annots - tp
            fn = label - tp
            for 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)/total
    mi = np.array(mi_list)/total
    smin = np.min(np.sqrt(ru * ru + mi * mi))
    return smin