理解MDAnalaysis DensityAnalysis中的密度对象

Question

我很难理解密度对象从DensityAnalysis类产生的确切结果。文档可以找到这里。

运行代码不是问题，而是了解密度对象到底产生了什么以及如何解释信息。

什么是"(113,113,113)回收箱“，我已经看到了MDAnalysis用户指南的例子，但我仍然不明白这是什么或如何解释它。

    from MDAnalysis.analysis.density import Density
    from MDAnalysis.analysis.density import DensityAnalysis
    from MDAnalysis import *
    import numpy as np

    PDB = '/Users/joveyosagie/Desktop/1vmd6lu7.pdb'
    DCD = '/Users/joveyosagie/Desktop/1vmd6lu7.dcd'
    u = Universe(PDB, DCD)
    protein = u.select_atoms('protein')
    OH2 = u.select_atoms('name OH2')

    OH2 = u.select_atoms('name OH2')    #select for water atoms
    D = DensityAnalysis(OH2, delta = 1.0) # each bin in histogram has size of 1 Angstrom
    D.run()

    D.density

[Out]Density density with (113, 113, 113) bins

Answer 1

MDAnalysis.analysis.density.Density对象保存用DensityAnalysis生成的三维直方图。密度产生的方式是通过计算粒子在一个小空间区域中出现的频率--体积元素或体素(带有长1的正交盒，但Density.delta精确地告诉你体素尺寸)或称为"bin“。我们对一个轨道的所有帧进行平均，并将计数归一化以得到一个密度(每体积粒子)。形状为NumPy的原始(num_bins_x, num_bins_y, num_bins_z)数组可作为Density.grid访问。

密度与原仿真的坐标系相关联。因此，我们还需要知道体素网格的来源(Density.origin保存这一信息)。使用origin、delta和数组的形状，我们现在可以计算出每个垃圾箱位于空间的位置。Density.edges属性提供沿x、y和z轴的bin边的值。例如，edges = [np.array(-2.5, -0.5, 1.5, 3.5]), np.array([0., 1., 2.]), np.array([2.5, 4.5])]属于具有形状(3, 2, 1)和delta = np.array(2.0, 1.0, 2.0])的网格。左下角的垃圾桶位于原点(-1.5, 0.5, 3.5) (原点在垃圾桶的中心)，包含坐标为-2.5≤x< -0.5、0≤y< 1和2.5≤z< 4.5的点。

该类包含更改密度存储单元(即Density.convert_density() )的方法。此方法通过将存储在grid中的值乘以适当的因子来更改底层数据。

其他方法是从构成gridData.core.Grid基础的Density类继承的。有关这些类的其他操作，请参阅GridDataFormats文档。例如，可以将两个Density对象作为numpy数组来处理，并对它们执行算法，例如减去它们以获得差异密度。

例子:比较水密度

如果密度是在相同的坐标系上生成的(即，具有相同的边)，则可以减去它们。

让我们比较两个模拟的水密度，看看配体对水的作用：

• apo (无配体)：u_apo宇宙
• holo (配体)：u_holo

首先，将轨迹叠加在一个共同的参考结构上，这样蛋白质就处于同一个坐标系中。您可以使用MDAnalysis.analysis.align.AlignTraj，如用AlignTraj对齐弹道中所述，也可以在用动态变换对中、对齐和使分子完整。下的“密度分析用户指南”中查看更详细的说明。

然后，我们需要确保我们的密度在同一个坐标系中被分析。

• 找到共同的参考中心。
• 在两种密度中设置相同数量的回收箱。我们假设一个30、x、30、x、30、多维数据集就足够了，但是您必须计算出正确的维度才能进行实际的分析。

protein_apo = u_apo.select("protein")
gridcenter = protein_apo.center_of_mass()  # should be the same as for holo

# select water oxygens for density
water_apo = u_apo.select_atoms("resname SOL and name OW")  # adjust for your simulation
water_holo = u_holo.select_atoms("resname SOL and name OW")

# perform the analysis
D_apo = DensityAnalysis(u_apo, gridcenter=gridcenter, xdim=30, ydim=30, zdim=30).run(verbose=True)
D_holo = DensityAnalysis(u_apo, gridcenter=gridcenter, xdim=30, ydim=30, zdim=30).run(verbose=True)

# work with the density instances
density_apo = D_apo.results.density
density_holo = D_holo.results.density

# convert to units in which water at standard conditions is 1
# (see Density.convert_units() docs for more choices)
density_apo.convert_units("water")
density_holo.convert_units("water")

# compare densities
delta_density = density_holo - density_apo

print("max difference", delta_density.grid.max())
print("min difference", delta_density.grid.min())

# export to DX file for visualization
delta_density.export("delta_holo_apo.dx")

更多的帮助？

如果你有更多的问题，请看看如何参与MDAnalysis社区 -我们有一个不和谐的服务器和邮件列表，人们(用户和开发人员)乐意帮助和讨论。