Berechnen Sie die Protein-Kontaktordnung in Python

Question

Protein contact order (CO) ist von Lokalität von Inter-Rückstandskontakte. Das CO korreliert auch mit der Faltungsgeschwindigkeit von Proteinen. Höhere Kontaktordnungen deuten auf längere Faltungszeiten hin, und eine geringe Kontaktordnung wurde als ein Prädiktor für eine mögliche Abwärtsfaltung oder Proteinfaltung, die ohne eine freie Energiebarriere auftritt, vorgeschlagen.

Es gibt einen Inline-Webserver und ein Perl-Skript zur Berechnung der CO here.

Gibt es eine Möglichkeit, das CO in Python zu berechnen?

Answer 1

Sie können die Tatsache nutzen, dass Sie wissen, dass Ihr Koordinatenfeld immer (N,3) in Form ist, so dass Sie eine Array-Erstellung und Aufruf np.dot loswerden können, die für wirklich kleine Arrays wie die hier weniger effizient ist. so tun, ich habe Ihre Funktion als abs_contact_order2 neu geschrieben:

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) 

Und dann Timings:

%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 

ich die print-Anweisungen aus den Funktionen übernommen haben, die Sie laufen zu lassen Numba Jit in nopython Modus. Wenn Sie diese Informationen wirklich benötigen, gebe ich alle erforderlichen Daten zurück und schreibe einen dünnen Wrapper, der die Rückgabewerte prüft und Debug-Informationen nach Bedarf ausgibt.

Beachten Sie auch, dass Sie beim Timing von Numba-Funktionen zuerst den Jit "aufwärmen" sollten, indem Sie die Funktion außerhalb der Timing-Schleife aufrufen. Wenn Sie die Funktion oft aufrufen, ist die Zeit, die Numba benötigt, um die Funktion zu verwenden, größer als die Zeit für einen einzelnen Anruf, so dass Sie nicht wirklich einen guten Hinweis darauf erhalten, wie schnell der Code ist. Wenn Sie jedoch die Funktion nur einmal aufrufen und die Startzeit wichtig ist, setzen Sie die Zeitmessung so fort, wie Sie sind.

Update:

Sie können dies auf ein wenig weiter beschleunigen, da Sie Schleifen über (i, j) und (j, i) Paare und sie sind symmetrisch:

@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) 

und timing:

%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

In der Tat gibt es! diese die ausgezeichnete MDTraj verwendet, ist das kein Problem:

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 mit:

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) 

Ohne numba diese 16s nimmt und nur durch eine einzige @jit Zugabe der Zeit zu ~ 1s reduziert wird.

Der ursprüngliche Perl-Skript dauert ca. 8s

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

Die Testdatei sowie ein jupyter Notebook here ist.