Рассчитать порядок контактов белка в Python

Question

Протеин contact order (CO) является локальным контактом между остатками. СО коррелирует также со скоростью сгибания белков. Более высокие порядки контактов указывают на более длительное время сгибания, и был предложен низкий порядок контактов в качестве предсказателя потенциального сгибания вниз или сгибание белка, которое происходит без барьера свободной энергии.

Для расчета CO here имеется встроенный веб-сервер и скрипт perl.

Есть ли способ вычислить CO в python?

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 

Я взял операторы печати из функций, которая позволяет запустить Numba jit в режиме nopython. Если вам действительно нужна эта информация, я бы вернул все необходимые данные и написал тонкую оболочку, которая проверяет возвращаемые значения и печатает любую отладочную информацию по мере необходимости.

Также обратите внимание, что при использовании функций Numba, вы должны «разогреть» jit сначала, вызывая функцию вне цикла синхронизации. Часто, если вы собираетесь вызывать функцию много раз, время, которое требуется Numba для jit функции, больше, чем время для одного вызова, поэтому вы действительно не получаете хороший показатель того, насколько быстро код. Если, однако, вы будете только один раз вызывать функцию, и время запуска важно, продолжайте синхронизировать ее так, как вы.

Update:

Вы можете ускорить это немного дальше, так как вы цикл над (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) 

Run с помощью:

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 время снижается до ~ 1 сек.

оригинальный сценарий Perl занимает около 8s

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

тест-файл, а как jupyter ноутбук here.