如何利用MDAnalysis计算蛋白质质量中心？

Question

我的处境有点不寻常。有七个不同的蛋白质存储在一个文件中，根据它们的残基名称。每种蛋白质都有不同的序列长度。现在，我需要计算每个蛋白质的质量中心，生成一个时间序列数据，我知道如何处理单个蛋白质，但不需要多个蛋白质系统。对于单一蛋白质，我可以这样做：

import MDAnalysis as mda
import numpy as np

u = mda.Universe('lp400start.gro')
u1 = mda.Merge(u.select_atoms("not resname W and not resname WF and not resname ION"))
u1.load_new('lp400.xtc')
protein = u1.select_atoms("protein")
arr = np.empty((protein.n_residues, u1.trajectory.n_frames, 3))

for ts in u.trajectory:
    arr[:, ts.frame] = protein.center_of_mass(compound='residues')

数据文件是可公开使用的这里。可以使用grep "^ *1[A-Z]" -B1 lp400final.gro | grep -v "^ *1[A-Z]"检查拓扑文件中的剩余序列，输出如下：

 38ALA     BB   74  52.901  33.911   6.318
--
   38ALA     BB  148  41.842  29.381   7.211
--
  137GLY     BB  455  36.756   4.287   3.284
--
  137GLY     BB  762  44.721  60.377   3.112
--
  252HIS    SC3 1368  28.682  37.936   6.727
--
  252HIS    SC3 1974  18.533  46.506   6.314
--
  576PHE    SC3 3263  48.937  38.538   4.013
--
  576PHE    SC3 4552  18.513  25.948   3.800
--
 1092PRO    SC1 6470  42.510  40.992   6.775
--
 1092PRO    SC1 8388  14.709   4.759   6.370
--
 1016LEU    SC110524  57.264  56.308   2.632
--
 1016LEU    SC112660  50.716  14.698   2.728
--
 1285LYS    SC215345   0.793  33.529   1.509

第一个蛋白质的序列长度为38个残基，它有自己的一个拷贝，然后是第二个蛋白质，等等。现在我想让每个蛋白质的COM在每个时间框架，并将它构建成一个时间序列。除了蛋白质之外，拓扑文件还包含DPPC粒子。Could someone help me how to do this?，谢谢！

为了确保输出轨迹是正确的，它看起来类似于这个在这里输入链接描述

Answer 1

我会从TPR文件中加载系统来维护债券信息。然后MDAnalysis可以确定片段(即你的蛋白质)。然后循环遍历片段以确定COM时间序列：

import MDAnalysis as mda
import numpy as np

# files from https://doi.org/10.5281/zenodo.846428
TPR = "lp400.tpr"
XTC = "lp400.xtc"

# build reduced universe to match XTC
# (ignore warnings that no coordinates are found for the TPR)
u0 = mda.Universe(TPR)
u = mda.Merge(u0.select_atoms("not resname W and not resname WF and not resname ION"))
u.load_new(XTC)

# segments (exclude the last one, which is DPPC and not protein)
protein_segments = u.segments[:-1]

# build the fragments
# (a dictionary with the key as the protein name -- I am using an
# OrderedDict so that the order is the same as in the TPR)
from collections import OrderedDict
protein_fragments = OrderedDict((seg.segid[6:], seg.atoms.fragments) for seg in protein_segments)

# analyze trajectory (with a nice progress bar)
timeseries = []
for ts in mda.log.ProgressBar(u.trajectory):
    coms = []
    for name, proteins in protein_fragments.items():
        # loop over all the different proteins;
        # unwrap to get the true COM under PBC (double check!!)
        coms.extend([p.center_of_mass(unwrap=True) for p in proteins]) 
    timeseries.append(coms)
timeseries = np.array(timeseries)

备注

• 再次检查unwrap=True是否在做正确的事情(而且这是必要的--我没有检查是否有任何蛋白质被跨周期分拆)。
• 展开是缓慢的；如果您不需要它，它将运行得更快。
• 得到的数组是一个具有(N_timesteps, M_proteins, 3)形状的三维数组，即(10001, 14, 3)。
• protein_fragments的含量为OrderedDict(<原子组与74 atoms>，<原子组与74 atoms>)，('OMPA‘)，(<原子组与307 atoms>，<原子组与307 atoms>)，('OMPG'，<原子组与606 atoms>，<原子组与606 atoms>)，('BTUB'，<原子组与1289 atoms>，<原子组与1289 atoms>))，('ATPS'，<原子组与1918年atoms>，)，('GLPF'，(，)，(‘atoms>’，，))