如何读取pdb文件并使用python执行基于FFT的对接？

Question

假设我有两个PDB文件(其中一个如下所示)

ATOM      1  N   MET A   1      66.104  56.583 -35.505  
ATOM      2  CA  MET A   1      66.953  57.259 -36.531  
ATOM      3  C   MET A   1      67.370  56.262 -37.627  
ATOM      4  O   MET A   1      67.105  55.079 -37.531  
ATOM      5  CB  MET A   1      68.227  57.852 -35.867  
ATOM      6  CG  MET A   1      67.848  58.995 -34.899  
ATOM      7  SD  MET A   1      66.880  58.593 -33.421  
....      .  ..  ... .   .      ......  ......  ......
....      .  ..  ... .   .      ......  ......  ......

可以使用以下脚本在python中读取此文件。

import sys
x=[];y=[];z=[]
res=[]
Nr=0
for fn in sys.argv[1:]:
   f=open(fn,'r')
   while 1:
      line=f.readline()
      if not line: break
      if line[0:6]=='ATOM  ' :
         rx=float(line[30:38]);ry=float(line[38:46]);rz=float(line[46:54])
         if line[21]=='A' :
            x.append(rx); y.append(ry); z.append(rz)
            Nr=Nr+1
            res.append(line[17:20])
   for i in range(1,Nr-1):
      print fn, i, res[i], x[i], y[i], z[i]
f.close

现在，我想生成N*N*N维数的网格，并在网格上旋转和平移分子。利用FFT (快速傅里叶变换)实现旋转和平移。

我试着写了如下的东西

import numpy as np
import fftw as fft

class Grid3D(object):

   def __init__(self, grid_dimension):
      x = y = z = grid_dimension
      self.grid = np.zeros([x, y, z], dtype=float)

所有这些实际上是用三维网格和FFT实现两个分子的对接。我想知道怎样做得更好还是更好？

Answer 1

回答您的第一个问题，“如何读取pdb文件”

如果您想最终得到一个numpy数组，您可以使用numpy.genfromtxt，它非常好，而且比循环读取更容易实现和使用。它对文件的间距也更健壮，等等。

import numpy as np
data = np.genfromtxt('filename.txt',
        names = 'ATOM,index,res,MET,A,count,x,y,z',
        dtype=['S4',int,'S2','S3','S1',int,float,float,float])

现在，data是一个numpy“结构化数组”，可以按以下方式轻松访问：

In [13]: data
Out[13]: 
array([('ATOM', 1, 'N', 'MET', 'A', 1, 66.104, 56.583, -35.505),
       ('ATOM', 2, 'CA', 'MET', 'A', 1, 66.953, 57.259, -36.531),
       ('ATOM', 3, 'C', 'MET', 'A', 1, 67.37, 56.262, -37.627),
       ('ATOM', 4, 'O', 'MET', 'A', 1, 67.105, 55.079, -37.531),
       ('ATOM', 5, 'CB', 'MET', 'A', 1, 68.227, 57.852, -35.867),
       ('ATOM', 6, 'CG', 'MET', 'A', 1, 67.848, 58.995, -34.899),
       ('ATOM', 7, 'SD', 'MET', 'A', 1, 66.88, 58.593, -33.421)], 
      dtype=[('ATOM', 'S4'), ('index', '<i8'), ('el', 'S2'), ('MET', 'S3'), ('A', 'S1'), ('count', '<i8'), ('x', '<f8'), ('y', '<f8'), ('z', '<f8')])

In [14]: data['x']
Out[14]: array([ 66.104,  66.953,  67.37 ,  67.105,  68.227,  67.848,  66.88 ])

In [15]: data['y']
Out[15]: array([ 56.583,  57.259,  56.262,  55.079,  57.852,  58.995,  58.593])

In [16]: data['index']
Out[16]: array([1, 2, 3, 4, 5, 6, 7])

In [17]: data[3]
Out[17]: ('ATOM', 4, 'O', 'MET', 'A', 1, 67.105, 55.079, -37.531)