对齐结构

载入文件

我们在此将要处理的测试文件是腺苷酸激酶(AdK)的运动轨迹，这是一种磷酸转移酶酶类。(BDPW09)这些轨迹展示了从闭合构象到开放构象的转变过程。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
import MDAnalysis as mda
from MDAnalysis.analysis import align
from MDAnalysis.tests.datafiles import CRD, PSF, DCD, DCD2
## import nglview as nv

import warnings
## suppress some MDAnalysis warnings about writing PDB files
warnings.filterwarnings('ignore')

adk_open = mda.Universe(CRD, DCD2)
adk_closed = mda.Universe(PSF, DCD)

## adk_open_view = nv.show_mdanalysis(adk_open)
## adk_open_view

## adk_closed_view = nv.show_mdanalysis(adk_closed)
## adk_closed_view

目前，这些蛋白质并未相互对齐。当将闭合构象与开放构象进行比较时，这种差异会更加明显。下面，我们将adk_open设为最后一个帧，并在合并后的宇宙中查看每个蛋白质的相对位置。

1
2
3
4
adk_open.trajectory[-1] ## last frame
merged = mda.Merge(adk_open.atoms, adk_closed.atoms)
## merged_view = nv.show_mdanalysis(merged)
## merged_view

使用align.alignto对齐结构

(API 文档)“alignto”函数会将移动的“AtomGroup”与目标“AtomGroup”进行对齐，以使粒子位置之间的“均方根偏差”(RMSD)最小化(有关 RMSD 的解释，请参阅链接的笔记本)。该函数返回(旧的 RMSD 值，新的 RMSD 值)。默认情况下(match_atoms=True)，它会尝试根据质量来匹配移动结构和参考结构之间的原子。

1
2
3
4
5
6
7
8
9
rmsds = align.alignto(adk_open,  ## mobile
                      adk_closed,  ## reference
                      select='name CA', ## selection to operate on
                      match_atoms=True) ## whether to match atoms
print(rmsds)
## (np.float64(21.712154435976014), 6.817293751703919)

## aligned_view = nv.show_mdanalysis(mda.Merge(adk_open.atoms, adk_closed.atoms))
## aligned_view

然而，您可能希望采用一种结构，其中粒子的质量没有明确的对应关系。例如，您可以将原子级蛋白质的α碳原子与粗粒化结构的主链节点对齐。下面，我们以一个稍微牵强的例子来说明，即将214个α碳原子与参考结构的前214个原子进行对齐。在这种情况下，我们需要将match_atoms设置为False，否则对齐操作将会出错。

1
2
3
4
5
6
7
8
9
10
rmsds = align.alignto(adk_open.select_atoms('name CA'),  ## mobile
                      adk_closed.atoms[:214],  ## reference
                      select='all', ## selection to operate on
                      match_atoms=False) ## whether to match atoms
print(rmsds)
## (np.float64(18.991465038265208), 16.603704620787127)

## shifted_aligned_view = nv.show_mdanalysis(mda.Merge(adk_open.atoms, adk_closed.atoms))
## shifted_aligned_view

当我们调整结构时，位置会暂时固定下来。如果我们将镜头切换到adk_open的第一帧，然后再切换回最后一帧，就会发现它已经回到了原来的位置。

1
2
3
4
5
6
adk_open.trajectory[0] ## set to first frame
adk_open.trajectory[-1] ## set to last frame
## < Timestep 101 >

## reset_view = nv.show_mdanalysis(mda.Merge(adk_open.atoms, adk_closed.atoms))
## reset_view

您可以将对齐后的位置信息写入一个 PDB 文件，并创建一个新的Universe模型，以此来保存这些位置信息。

1
2
3
4
5
align.alignto(adk_open, adk_closed, select='name CA')
adk_open.atoms.write('aligned.pdb')

## from_file_view = nv.show_mdanalysis(mda.Universe('aligned.pdb'))
## from_file_view

对齐轨迹

载入文件

我们将在这里使用的测试文件是腺苷酸激酶 (AdK) 的轨迹，这是一种磷酸转移酶。(BDPW09)轨迹对从封闭构象到开放构象的转变进行采样。

1
2
3
4
5
6
7
8
9
10
11
import MDAnalysis as mda
from MDAnalysis.analysis import align
from MDAnalysis.tests.datafiles import CRD, PSF, DCD, DCD2
## import nglview as nv

import warnings
## suppress some MDAnalysis warnings when writing PDB files
warnings.filterwarnings('ignore')

adk_open = mda.Universe(CRD, DCD2)
adk_closed = mda.Universe(PSF, DCD)

使用AlignTraj对齐轨迹

align.alignto用于对结构或轨迹的某一帧进行对齐，而align.AlignTraj(API文档)则能高效地将整个轨迹与参考值进行对齐。与大多数其他分析模块不同，AlignTraj允许您将分析结果输出到文件中。这是因为当通过从文件加载的方式来创建宇宙时，对于逐帧(动态)信息所做的更改在帧发生改变时不会保持不变。如果轨迹未被写入文件，或者未被加载到内存中(如下所示)，那么AlignTraj将不会产生任何效果。

1
2
3
4
5
6
7
8
9
10
11
align.AlignTraj(adk_closed,  ## trajectory to align
                adk_open,  ## reference
                select='name CA',  ## selection of atoms to align
                filename='aligned.dcd',  ## file to write the trajectory to
                match_atoms=True,  ## whether to match atoms based on mass
               ).run()
## merge adk_closed and adk_open into the same universe
merged1 = mda.Merge(adk_closed.atoms, adk_open.atoms)

## merged1_view = nv.show_mdanalysis(merged1)
## merged1_view

如您所见，adk_closed和adk_open这两条轨迹看起来仍然是一样的。然而，当我们从aligned.dcd文件中加载已对齐的轨迹时，就会发现它们是重叠在一起的：

1
2
3
4
5
6
7
aligned = mda.Universe(PSF, 'aligned.dcd')
aligned.segments.segids = ['Aligned']  ## rename our segments
adk_open.segments.segids = ['Open']  ## so they're coloured differently
merged2 = mda.Merge(aligned.atoms, adk_open.atoms)

## merged2_view = nv.show_mdanalysis(merged2)
## merged2_view

如果不想写文件，也可以选择将整个轨迹加载到内存中。(这并不总是可行的，具体取决于您的轨迹有多大，以及您的设备有多少内存，在这种情况下，如上所述将对齐轨迹写入文件会更有效)。您可以通过以下两种方式之一完成此作：

  1. 将整个轨迹加载到内存中，
  2. 在AlignTraj中使用in_memory=True。

复制坐标到新的Universe

MDAnalysis.Merge 操作不会自动加载轨迹中的坐标数据。我们可以自己进行这项操作。下面，我们将复制对齐后的体系中的 98 个帧的坐标信息。

1
2
3
4
5
6
7
8
9
10
11
12
from MDAnalysis.analysis.base import AnalysisFromFunction
import numpy as np
from MDAnalysis.coordinates.memory import MemoryReader

def copy_coords(ag):
    return ag.positions.copy()

aligned_coords = AnalysisFromFunction(copy_coords,
                                      aligned.atoms).run().results

print(aligned_coords['timeseries'].shape)
## (98, 3341, 3)

相比之下，我们将让开放构象保持不变。

1
2
3
adk_coords = adk_open.atoms.positions.copy()
adk_coords.shape
## (3341, 3)

因为对齐后的Universe共有 98 个帧，所以我们复制了 adk_open 位置的坐标，并将它们进行堆叠。

1
2
3
adk_traj_coords = np.stack([adk_coords] * 98)
adk_traj_coords.shape
(98, 3341, 3)

我们通过使用 np.hstack 方法将 aligned_coords 和 adk_traj_coords 在第二个轴上进行堆叠连接，并将这些坐标加载到内存中的 merged2 宇宙中。

1
2
3
4
5
6
7
merged_coords = np.hstack([aligned_coords['timeseries'],
                           adk_traj_coords])
merged2.load_new(merged_coords, format=MemoryReader)
## <Universe with 6682 atoms>

## m2_view = nv.show_mdanalysis(merged2)
## m2_view

在线notebook无法展示分子的运动轨迹，但在这里您可以使用 nglview.contrib.movie.MovieMaker 来制作该轨迹的GIF图片。这需要您先安装 moviepy 插件。

1
2
3
4
5
6
7
8
## from nglview.contrib.movie import MovieMaker
## movie = MovieMaker(
##     m2_view,
##     step=4,  ## only render every 4th frame
##     output='merged.gif',
##     render_params={"factor": 3},  ## set to 4 for higher quality
## )
## movie.make()

将轨迹写入文件

最后，我们还可以将这个将新轨迹保存一个文件中。

1
2
3
with mda.Writer('aligned.xyz', merged2.atoms.n_atoms) as w:
    for ts in merged2.trajectory:
        w.write(merged2.atoms)

对齐轨迹本身

载入文件

我们在此将要处理的测试文件是腺苷酸激酶(AdK)的运动轨迹，这是一种磷酸转移酶酶类。(BDPW09)这些轨迹展示了从闭合构象到开放构象的转变过程。

1
2
3
4
5
6
import MDAnalysis as mda
from MDAnalysis.analysis import align, rms
from MDAnalysis.tests.datafiles import PSF, DCD

mobile = mda.Universe(PSF, DCD)
ref = mda.Universe(PSF, DCD)

对齐轨迹到第一帧

align.alignto 用于对结构或轨迹的框架进行对齐，而 align.AlignTraj(API 文档)则能高效地将整个轨迹与参考轨迹进行对齐。
我们首先会检查未对齐轨迹的均方根偏差(RMSD)值，以便能够进行结果比较(请查看链接中的笔记本以了解RMSD的解释)。下面的代码通过索引最后一个时间步长将移动轨迹设置为最后一个帧，通过索引第一个时间步长将ref设置为第一个帧，并计算α-碳位置之间的均方根偏差。

1
2
3
4
5
6
7
8
mobile.trajectory[-1]  ## set mobile trajectory to last frame
ref.trajectory[0]  ## set reference trajectory to first frame

mobile_ca = mobile.select_atoms('name CA')
ref_ca = ref.select_atoms('name CA')
unaligned_rmsd = rms.rmsd(mobile_ca.positions, ref_ca.positions, superposition=False)
print(f"Unaligned RMSD: {unaligned_rmsd:.2f}")
## Unaligned RMSD: 6.84

现在我们可以对轨迹进行对齐了。我们已经将ref设置为第一帧。在下面的单元格中，我们将轨迹的位置加载到内存中，以便在Python中对其进行修改。

1
aligner = align.AlignTraj(mobile, ref, select='name CA', in_memory=True).run()

如果您没有足够的内存来完成此操作，请将轨迹信息写入文件，然后将其重新加载到 MDAnalysis 中（取消注释下面的单元格）。否则，您无需运行该操作。

1
2
3
## aligner = align.AlignTraj(mobile, ref, select='backbone',
##                           filename='aligned_to_first_frame.dcd').run()
## mobile = mda.Universe(PSF, 'aligned_to_first_frame.dcd')

现在我们可以看到，均方根偏差已经有所降低（幅度较小）。

1
2
3
4
5
6
7
8
9
mobile.trajectory[-1]  ## set mobile trajectory to last frame
ref.trajectory[0]  ## set reference trajectory to first frame

mobile_ca = mobile.select_atoms('name CA')
ref_ca = ref.select_atoms('name CA')
aligned_rmsd = rms.rmsd(mobile_ca.positions, ref_ca.positions, superposition=False)

print(f"Aligned RMSD: {aligned_rmsd:.2f}")
## Aligned RMSD: 6.81

对齐轨迹到第三帧

我们可以将轨迹调整到任何帧中：比如第三帧。这个过程大致相同，只是我们需要通过索引第三个时间步来将ref设置为第三帧。

1
2
3
4
5
6
7
mobile.trajectory[-1]  ## set mobile trajectory to last frame
ref.trajectory[2]  ## set reference trajectory to third frame

aligned_rmsd_3 = rms.rmsd(mobile.atoms.positions, ref.atoms.positions, superposition=False)

print(f"Aligned RMSD: {aligned_rmsd_3:.2f}")
## Aligned RMSD: 6.73

计算原子结构的均方根偏差值

翻译https://userguide.mdanalysis.org/stable/examples/analysis/alignment_and_rms/aligning_structure_to_another.html