0

假设我有两个 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)

这一切实际上是为了使用 3d 网格和 FFT 对两个分子进行对接。我想知道如何进一步或更好的方法？

4

1 回答 1

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)
于 2013-04-01T15:53:25.603 回答