用Python计算蛋白质接触顺序

Question

蛋白质接触顺序 ( Protein，CO)是残基间接触的局部性.CO也与蛋白质的折叠速度有关。较高的接触级数表明折叠时间较长，而低接触序被认为是潜在的下坡折叠或无自由能垒发生的蛋白质折叠的预测因子。

有一个内联the服务器和一个perl脚本来计算CO 这里。

有没有方法计算蟒蛇的CO？

Answer 1

您可以利用这样一个事实，即您知道您的坐标数组总是形状为(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

我已经从函数中删除了print语句，它允许您在nopython模式下运行Numba。如果您确实需要这些信息，我将返回所有必要的数据，编写一个细包装器，检查返回值并根据需要打印任何调试信息。

还请注意，当计时Numba函数时，您应该首先通过调用定时循环之外的函数来“热身”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

Answer 2

的确有！使用优秀的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，这需要16s，只要添加一个@jit，时间就会减少到~1s。

最初的perl脚本大约需要8s。

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

测试文件和jupyter笔记本都是这里。