计算 Python 中的蛋白质接触顺序

Question

蛋白质contact order (CO) 是残基间接触的局部。 CO 还与蛋白质的折叠速度相关。较高的接触阶数表示较长的折叠时间，低接触阶数已被建议作为潜在下坡折叠或在没有自由能垒的情况下发生的蛋白质折叠的预测指标。

有一个内联网络服务器和一个 perl 脚本来计算 CO here。

有没有办法计算 python 中的二氧化碳？

Answer 1

确实有！使用优秀的MDTraj这不是问题:

import numpy as np
import mdtraj as md

@jit
def abs_contact_order(xyz, atoms_residue, cutoff_nm=.6):
    """Return the absolute contact order."""
    contact_count = 0
    seq_distance_sum = 0
    cutoff_2 = cutoff_nm*cutoff_nm
    N = len(atoms_residue)
    for i in xrange(N):
        for j in xrange(N):
            seq_dist = atoms_residue[j] - atoms_residue[i]
            if seq_dist > 0:
                d = xyz[j] - xyz[i]
                if np.dot(d, d) < cutoff_2:
                    seq_distance_sum += seq_dist 
                    contact_count += 1


    if contact_count==0.:
        print("Warning, no contacts found!")
        return 0.

    print("contact_count: ", contact_count)
    print("seq_distance_sum: ", seq_distance_sum)
    return seq_distance_sum/float(contact_count)

运行使用：

traj = md.load("test.pdb")
print("Atoms: ", traj.n_atoms)
seq_atoms = np.array([a.residue.resSeq for a in traj.top.atoms], dtype=int)

abs_co = abs_contact_order(traj.xyz[0], seq_atoms, cutoff_nm=0.60)

print("Absolute Contact Order: ", abs_co)
print("Relative Contact Order: ", abs_co/traj.n_residues)

如果没有 numba，这需要 16 秒，只需添加一个 @jit，时间就会减少到 ~1 秒。

原来的perl脚本大约需要8s

perl contactOrder.pl -r -a -c 6 test.pdb

测试文件和 jupyter 笔记本是 here。

Answer 2

您可以利用您知道坐标数组的形状始终为 (N,3) 的事实，因此您可以摆脱数组创建和调用 np.dot，这会降低效率对于像这里的那些非常小的阵列。这样做，我将您的函数重写为 abs_contact_order2:

from __future__ import print_function, division

import mdtraj as md
import numpy as np
import numba as nb

@nb.njit
def abs_contact_order(xyz, atoms_residue, cutoff_nm=.6):
    """Return the absolute contact order."""
    contact_count = 0
    seq_distance_sum = 0
    cutoff_2 = cutoff_nm*cutoff_nm
    N = len(atoms_residue)
    for i in xrange(N):
        for j in xrange(N):
            seq_dist = atoms_residue[j] - atoms_residue[i]
            if seq_dist > 0:
                d = xyz[j] - xyz[i]
                if np.dot(d, d) < cutoff_2:
                    seq_distance_sum += seq_dist 
                    contact_count += 1


    if contact_count==0.:
        return 0.

    return seq_distance_sum/float(contact_count)

@nb.njit
def abs_contact_order2(xyz, atoms_residue, cutoff_nm=.6):
    """Return the absolute contact order."""
    contact_count = 0
    seq_distance_sum = 0
    cutoff_2 = cutoff_nm*cutoff_nm
    N = len(atoms_residue)
    for i in xrange(N):
        for j in xrange(N):
            seq_dist = atoms_residue[j] - atoms_residue[i]
            if seq_dist > 0:
                d = 0.0
                for k in xrange(3):
                    d += (xyz[j,k] - xyz[i,k])**2
                if d < cutoff_2:
                    seq_distance_sum += seq_dist 
                    contact_count += 1


    if contact_count==0.:
        return 0.

    return seq_distance_sum/float(contact_count)

然后时间：

%timeit abs_co = abs_contact_order(traj.xyz[0], seq_atoms, cutoff_nm=0.60)
1 loop, best of 3: 723 ms per loop

%timeit abs_co = abs_contact_order2(traj.xyz[0], seq_atoms, cutoff_nm=0.60)
10 loops, best of 3: 28.2 ms per loop

我已经从函数中删除了打印语句，这样您就可以运行 nopython 模式下的 Numba jit。如果您真的需要这些信息，我会 return 所有必要的数据并编写一个薄包装程序来检查 return 值并根据需要打印任何调试信息。

另请注意，在为 Numba 函数计时时，您应该 "warm-up" 通过在计时循环之外调用该函数来首先执行 jit。通常，如果您要多次调用该函数，Numba jit 函数所花费的时间比单次调用所花费的时间要长，因此您无法真正了解代码的速度有多快。但是，如果您只打算调用该函数一次并且启动时间很重要，请继续按原样计时。

更新：

您可以进一步加快速度，因为您正在遍历 (i,j) 和 (j,i) 对并且它们是对称的：

@nb.njit
def abs_contact_order3(xyz, atoms_residue, cutoff_nm=.6):
    """Return the absolute contact order."""
    contact_count = 0
    seq_distance_sum = 0
    cutoff_2 = cutoff_nm*cutoff_nm
    N = len(atoms_residue)
    for i in xrange(N):
        for j in xrange(i+1,N):
            seq_dist = atoms_residue[j] - atoms_residue[i]
            if seq_dist > 0:
                d = 0.0
                for k in xrange(3):
                    d += (xyz[j,k] - xyz[i,k])**2
                if d < cutoff_2:
                    seq_distance_sum += 2*seq_dist 
                    contact_count += 2


    if contact_count==0.:
        return 0.

    return seq_distance_sum/float(contact_count)

和时间：

%timeit abs_co = abs_contact_order3(traj.xyz[0], seq_atoms, cutoff_nm=0.60)
100 loops, best of 3: 15.7 ms per loop

计算 Python 中的蛋白质接触顺序

Calculate protein contact order in Python

python

numba