Python通过坐标系变换实现分子结构的平移和旋转

背景

在最近做的课题中需要将蛋白质靶点的圆柱形空腔展开，由于分子的坐标本身不是沿坐标轴展开的，因此处理起来异常的麻烦。以我浅薄的知识判断，把蛋白-配体形成的腔体的主轴放到x轴上，在在y-z方向上做极坐标展开应该是一个简单可行的方法。

在实现上，把主轴平移到x轴和直接把主轴作为x轴是等效的。然而这里面最shaib的地方在于，网上的很多教程写的乱七八糟，对数学不好的我来说那真是吃屎一样难受，因此自己简单做了做整理了一下。

坐标变换原理

这里主要参考这一篇博客，注意这里的逆矩阵计算有些问题
图形学1-三维坐标系间的变换矩阵推导.
简单点说，变换前后的坐标满足这些关系式：

其中（x0，y0，z0）是新坐标系的原点在原坐标系的坐标，（x，y，z）是要变换的点在原坐标系的坐标。其他的符号可以参考原博客。这些关系式建立了这样的联系：新坐标系下的坐标加上坐标原点的平移量（比如x0）等于原坐标。这样就可以得到：

这里我们暂时把这个4×4的矩阵称为T，显然原坐标和T^-1的矢量积就是新坐标了。
因为numpy中的np.linalg.inv()已经帮我们解决了逆矩阵的问题，这里我们就不做推导了。

实现

体系的初始结构如上图所示。这里面首先要处理的问题是新坐标系的基向量的方向。这里我选择了1493和1495两个残基上的α碳的方向作为x轴，因为这两几个碳原子的结构比较稳定。随后选择1494号残基的α碳，把1494-αC和x轴的垂足的方向做y轴，把x轴、y轴的法向量做z轴。
垂足的计算可以参考这篇博客：三维空间：点到直线垂足坐标公式推导.
在确定了新的坐标轴后，计算矩阵并求逆，并将获得的矩阵和原坐标相乘即可。可以看到配体已经被移动到了新的坐标系的原点附近。

代码

整个流程的代码如下

#载入gro文件，这里我去除了水和离子，保存为dat格式
ResID,ResName,AtomName,x,y,z,OUTPUT = __readfile("sample.dat")

#找到1493和1495的α碳
a = []
for i in range(0,len(ResID)):
    if ResID[i] == 1493 and AtomName[i] == "CA" or ResID[i] == 1495 and AtomName[i] == "CA" :
        #print(ResName[i])
        a.append(i)
        #找到1494的α碳
    elif ResID[i] == 1494 and AtomName[i] == "CA":
        b = i

#把1493α碳记为x1，1495α碳记为x2，1494α碳记为x0
x1 = np.array((x[a[0]],y[a[0]],z[a[0]]))
x2 = np.array((x[a[1]],y[a[1]],z[a[1]]))
x0 = np.array((x[b],y[b],z[b]))
#获得第一个基向量u
u = x1-x2
#去取1493和1495的中点为新的坐标原点
P0 = ((x[a[0]]+x[a[1]])/2,(y[a[0]]+y[a[1]])/2,(z[a[0]]+z[a[1]])/2)

#计算垂足，公式太长分子分母分开写
fz = (x[a[1]]-x[b])*(x[a[0]]-x[a[1]])+ (y[a[1]]-y[b]) * (y[a[0]]-y[a[1]]) + (z[a[1]]-z[b])*(z[a[0]]-z[a[1]])
fm = (x[a[0]]-x[a[1]])*(x[a[0]]-x[a[1]])+ (y[a[0]]-y[a[1]])*(y[a[0]]-y[a[1]])+ (z[a[0]]-z[a[1]])*(z[a[0]]-z[a[1]])
k = -fz/fm
#Pn为垂足的坐标
Pn = ( k*(x[a[0]]-x[a[1]]) + x[a[1]] , k*(y[a[0]]-y[a[1]]) + y[a[1]] , k*(z[a[0]]-z[a[1]]) + z[a[1]] )
#垂足和1494碳的方向确定为新y轴，获得第二个基向量v
v = np.array(Pn) - np.array(x0)

#第三个基向量由u、v确定的平面计算得出
w = np.array([(u[1]*v[2]-v[1]*u[2]),(v[0]*u[2]-u[0]*v[2]),(u[0]*v[1]-v[0]*u[1])])

#基向量的长度等于1
u = u/np.linalg.norm(u)
v = v/np.linalg.norm(v)
w = w/np.linalg.norm(w)


#获得矩阵T
T = np.array([[u[0],u[1],u[2],P0[0]],
      [v[0],v[1],v[2],P0[1]],
      [w[0],w[1],w[2],P0[2]],
      [0,0,0,1]])

#获得逆矩阵
T_1 = np.linalg.inv(T)

#输出到gro文件
ff = open("test.gro","w")
print("Test File",file=ff)
print(" "+str(len(ResID)),file=ff)
for i in range(0,len(ResID)):
    coord0 = np.array([x[i],y[i],z[i],1])
    newcoord = np.dot(GG,coord0)#这里的计算的矢量积即为新坐标
    p1 = '%.3f' % newcoord[0]#gro文件取三位浮点数
    p2 = '%.3f' % newcoord[1]
    p3 = '%.3f' % newcoord[2]
    print(OUTPUT[i]+p1.rjust(8)+p2.rjust(8)+p3.rjust(8),file=ff)
print("  14.80330  11.16660  15.71149",file=ff)
ff.close()

对于只考虑坐标变换的同学，直接使用这个函数即可：

#输入：新坐标系的基向量u、v、w和新坐标原点P0
def TransfMatrix(u,v,w,P0):
    u = u/np.linalg.norm(u)
    v = v/np.linalg.norm(v)
    w = w/np.linalg.norm(w)
    T = np.array([[u[0],u[1],u[2],P0[0]],
          [v[0],v[1],v[2],P0[1]],
          [w[0],w[1],w[2],P0[2]],
          [0,0,0,1]])
    return np.linalg.inv(T)
#输出坐标系变换矩阵
#注意顺序！！
NewCoord = np.dot(T_1,OriginCoord)