Tanimoto系数计算方法优化的思路
Kotori Y 27 Posts
 2019-01-08 11:21 2022-08-20 23:16      

将原本的嵌套循环简化成一个循环，大幅缩短计算时间，提高效率~

作为本Blog的第一篇正式文章，就写写前一阵子，我优化一个嵌套for循环的思路吧。な、始まるよ

相关背景

雅卡尔指数

雅卡尔指数（Jaccard index），又称为并交比（Intersection over Union）、雅卡尔相似系数（Jaccard similarity coefficient），是用于比较样本集的相似性与多样性的统计量。雅卡尔系数能够量度有限样本集合的相似度，其定义为两个集合交集大小与并集大小之间的比例：

Intersection

Union

若样本A与样本B完全重叠，显然 $J(A,B) = 1$，所以有：$0 \leq J(A,B) \leq1$

Tanimoto 系数

Tanimoto 相似性（Tanimoto similarity），雅卡尔相似性经常被误解成Tanimoto 相似性，又或者说两者的界限其实很模糊？（好吧，其实我也弄不懂……)

Tanimoto 相似性$T_s(X,Y) = \frac{\sum_i\left({X_i}\land{Y_i}\right)}{\sum_i\left({X_i}\lor{Y_i}\right)}$   Tanimoto 距离$T_d = -\log_2 \left(T_s\left(X,Y\right)\right)$

如果通过bit vector比较雅卡尔相似性Tanimoto 相似性，那它们可以表示成：$f(A,B) = \frac{A \cdot B}{\vert A \vert^2 + \vert B \vert^2 - A \cdot B}$

其中，$A \cdot B = \sum_i A_i B_i = \sum_i\left(A_i \land B_i \right); \vert A \vert^2 = \sum_i A_i^2$

我们要计算的就是这个系数~

CATS 描述符

CATS(Chemically Advanced Template Search) 描述符是一种基于药效团的描述符，下面介绍一种计算方法：

我们定义五种原子类型：氢键供体（hydrogen-bond donor, D）、氢键受体（hydrogen-bond acceptor, A）、正电荷（positively charged, P）、负电荷（negatively charged, N）和亲脂性（lipophilic, L）。

Example

将二维平面的分子（A）转换成molecular graph(B)

那么一共可以产生15对原子类型对（DD, DA, DP, DN, DL, AA, AP, AN, AL, PP, PN, PL, NN, NL, LL），即$C_5^2 + 5$。对于每个被划分到上述五种类型的原子，我们从n = 1遍历到n = 10，统计以某个原子为中心，统计n个键的范围内15种原子类型对的数量。这样我们会的到10个频数分布直方图，用每个直方图的高度除以分子中非氢原子的数量，最后得到一个150维（$15 \times 10$）的特征向量，用以表征一个分子，供计算机识别。

MOE(2018.01)插件中计算得到的是一个210维的CATS 描述符，所以我猜测（ん？）MOE中定义了6种原子类型（$\left(C_6^2 + 6 \right) \times 10 = 210$)。而这次项目要计算1660个分子间Tanimoto相似性，即处理一个1660行，210列的二维数组。

实现思路

误区

MOE计算的CATS描述符（部分）如下：

CATS

可见，CATS中是含有0的，由此我们便陷入了一个误区——对于计算A, B两个分子的Tanimoto系数，我们要找出A, B两分子的CATS中都不为0特征，即找出A, B的交集以计算$A \cdot B$、$\vert A \vert^2$、$\vert B \vert^2$。为此，我想了很多方法，但效果甚微（其中有个要运行20+h……）。

解决方法

其实，我们根本不需要找出A, B两分子的CATS中都不为0特征。对于两个分子的特征向量，我们可以直接进行数组间的运算，即$arr1 \times arr2$， $arr1^2$，$arr2^2$，因为0在乘法运算中就是在做交集运算!

接下来，我们只需要解决排列组合问题即可。我在网上发现了product这个函数：

1
2
3
4
from itertools import product

for i in product([1,2,3],repeat= 2):
    print(i)

点击查看运行结果 
1
2
3
4
5
6
7
8
9
(1, 1)
(1, 2)
(1, 3)
(2, 1)
(2, 2)
(2, 3)
(3, 1)
(3, 2)
(3, 3)

虽然，乍一看重复计算了将近一倍（(1,2), (2,1)…），但它是按“顺序”经行遍历的，最后的结果刚好能输出成一个$n \times n$的矩阵。代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
from itertools import product
import pandas as pd
import numpy as np

def CalTanimoto(data):
    arr = np.array(data) #data为1660个分子的CATS描述符，dataframe型
    
    start = time.clock()
    
    values = [] #用于存放每个得到的Tanimoto相似度
    for pair in product(list(range(1660)), repeat=2): #进行排列组合运
            up = arr[pair[0]]*arr[pair[-1]] #获取两个特征向量的数组，进行相关运算，下同
            down = arr[pair[0]]**2 + arr[pair[-1]]**2
            value = up.sum()/(down.sum()-up.sum())#计算	
            values.append(value)
            
    Tani = np.array(values).reshape((1660,1660))
    
    end = time.clock()
    print(start-end)
	
    return Tani

总计用时：约240s

解决方法2

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
import numpy as np
import warnings
warnings.filterwarnings('ignore')

# 更新时间：11/09/20


def getSimilarity(data, binary=True):
    """
    
    Parameters
    ----------
    data : numpy.ndarray
        fingrtprint or descriptor in 2d-arrry.
    Yields
    ------
    value : numpy.ndarray
        the similarity between a molecule and the lib.
    """
    for arr in data:
        up = arr*data #获取两个特征向量的数组，进行相关运算，下同
        
        if binary:
            down = arr + data
        else:
            down = arr**2 + data**2
            
        value = up.sum(axis=1)/(down.sum(axis=1)-up.sum(axis=1))#计算公式
        yield value
   

def main(data, binary=True):
    """
    Parameters
    ----------
    data : numpy.ndarray
        DESCRIPTION.
    Returns
    -------
    Similarity matrix.
    """
    simi = getSimilarity(data, binary)
    try:
        simi = np.vstack(simi)
    except:
        simi = np.array(list(simi))
    
    return simi
        
        

if "__main__" == __name__:
    import pandas as pd
    data = np.random.randint(0,2,(1660,1660)) #fake data
    simi = main(data)

解决方法3

1
2
3
4
5
6
7
8
9
10
11
12
13
14
import numpy as np
from numpy.typing import ArrayLike


def tanimoto_3(array: ArrayLike) -> float:
    m = array.shape[0]

    mat = np.mat(array)
    mat2 = np.mat(array**2).sum(axis=1).repeat(m).reshape(m, -1)

    up = mat * mat.T
    down = mat2 + mat2.T

    return up.__array__() / (down - up).__array__()  # 原则上矩阵无法做分母，但numpy可以这么做

测试

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
def test0():  # 测试方法1
    np.random.seed(777)
    array = np.random.randint(0, 2, (1000, 1660))
    CalTanimoto(array)


def test1():  # 测试方法2
    np.random.seed(777)
    array = np.random.randint(0, 2, (1000, 1660))
    main(array)


def test2():  # 测试方法3
    np.random.seed(777)
    array = np.random.randint(0, 2, (1000, 1660))
    tanimoto_3(array)


if "__main__" == __name__:
    from timeit import timeit

    t0 = timeit("test0()", "from __main__ import test0", number=10)  # 82.92235608300001
    t1 = timeit("test1()", "from __main__ import test1", number=10)  # 38.50827858299999
    t2 = timeit("test2()", "from __main__ import test2", number=10)  # 6.618927875000011

相关链接

1.雅卡尔指数wiki

2.Jaccard index wiki

参考文献

  1. Reutlinger, M., Koch, C. P., Reker, D., Todoroff, N., Schneider, P., Rodrigues, T., & Schneider, G. (2013). Chemically advanced template search (CATS) for scaffold‐hopping and prospective target prediction for ‘orphan’molecules. Molecular informatics, 32(2), 133-138.
  2. Schneider, G., Neidhart, W., Giller, T., & Schmid, G. (1999). “Scaffold‐hopping” by topological pharmacophore search: a contribution to virtual screening. Angewandte Chemie International Edition, 38(19), 2894-2896.

后记

这篇文章到就差不多了，作为第一篇文章，写了好几天（才不是因为我懒呢……)。为了写这篇文章，也看了不少东西：MathJax语法，HTML，还有两篇文献。还参照了dalao的写作格式。另外，感觉能有更好的方法来替代排列组合运算，进一步提高效率。但是，等我想出改进方法，恐怕我都算完20+个矩阵了。最后，我想说：“挖坑一时爽，一直挖坑一直爽”。

后记2

将近两年后再看这篇博客，真是感慨万千，当初在老哥的带领下我励志要做技术博客，没有想到却用来记录琐碎无聊的日常，而带我走进前端开发的老哥已不在CBDD。不过话说回来，那时的我也算初露锋芒，写出了性能还凑合的解决方案，而那个“更好的方法”也在一年半后实现了，也算有始有终吧。

后记3

在原文章的发布后的44个月，我终于去掉了函数中不优雅的for，从最开始的枚举到现在小o所说的“矩阵的思维”。希望Tanimoto这个函数的性能还能进一步的提升。

  • Post title：Tanimoto系数计算方法优化的思路
  • Post author：Kotori Y
  • Create time：2019-01-08 11:21
  • Update time：2022-08-20 23:16
  • Post link：https://blog.iamkotori.com/2019/01/08/Tanimoto系数计算方法优化/
  • Copyright Notice：All articles in this blog are licensed under BY-NC-SA unless stating additionally.
 Comments