python 为什么这个源代码给出了不正确的RMSD值？

nfeuvbwi 于 2023-05-05 发布在 Python

关注(0)|答案(1)|浏览(180)

将PDB结构与Biopython对齐

下面的源代码是从上面的链接获得的：

import Bio.PDB    
    # Select what residues numbers you wish to align
    # and put them in a list
    start_id = 1
    end_id   = 70
    atoms_to_be_aligned = range(start_id, end_id + 1)
    
    # Start the parser
    pdb_parser = Bio.PDB.PDBParser(QUIET = True)
    
    # Get the structures
    ref_structure = pdb_parser.get_structure("reference", "c:/home/bbq_input_pdb/pdb1a6j.pdb")
    sample_structure = pdb_parser.get_structure("sample", "c:/home/bbq_output_pdb/1a6j.pdb")
    
    # Use the first model in the pdb-files for alignment
    # Change the number 0 if you want to align to another structure
    ref_model    = ref_structure[0]
    sample_model = sample_structure[0]
    
    # Make a list of the atoms (in the structures) you wish to align.
    # In this case we use CA atoms whose index is in the specified range
    ref_atoms = []
    sample_atoms = []
    
    # Iterate of all chains in the model in order to find all residues
    for ref_chain in ref_model:
      # Iterate of all residues in each model in order to find proper atoms
      for ref_res in ref_chain:
        # Check if residue number ( .get_id() ) is in the list
        if ref_res.get_id()[1] in atoms_to_be_aligned:
          # Append CA atom to list
          ref_atoms.append(ref_res['CA'])
    
    # Do the same for the sample structure
    for sample_chain in sample_model:
      for sample_res in sample_chain:
        if sample_res.get_id()[1] in atoms_to_be_aligned:
          sample_atoms.append(sample_res['CA'])
    
    # Now we initiate the superimposer:
    super_imposer = Bio.PDB.Superimposer()
    super_imposer.set_atoms(ref_atoms, sample_atoms)
    super_imposer.apply(sample_model.get_atoms())
    
    # Print RMSD:
    print("RMSD = ", super_imposer.rms)
    
    ## Save the aligned version of 1UBQ.pdb
    #io = Bio.PDB.PDBIO()
    #io.set_structure(sample_structure)
    #io.save("1UBQ_aligned.pdb")

此源代码比较两个PDB蛋白质文件并计算它们的RMSD值。
然而，我观察到：
1.这个工具给出的值与PyMOL给出的值不同
1.即使我们将一个文件与其自身进行比较，它也会显示值-在PyMOL中，不可能将一个文件与其自身进行比较
谁能告诉我bug在哪里，因为我找不到它？

python

来源：https://stackoverflow.com/questions/74152346/why-is-this-source-code-giving-incorrect-rmsd-value

1条答案

按热度按时间

njthzxwz1#

我不是这里的Maven，我只是想从代码中找出答案;参见https://github.com/biopython/biopython/blob/master/Bio/PDB/Superimposer.py#L43：

import numpy as np

from Bio.SVDSuperimposer import SVDSuperimposer -------------------------> HERE !!!!!!!!
from Bio.PDB.PDBExceptions import PDBException

class Superimposer:
    """Rotate/translate one set of atoms on top of another to minimize RMSD."""

    def __init__(self):
        """Initialize the class."""
        self.rotran = None
        self.rms = None

    def set_atoms(self, fixed, moving):
        """Prepare translation/rotation to minimize RMSD between atoms.
        Put (translate/rotate) the atoms in fixed on the atoms in
        moving, in such a way that the RMSD is minimized.
        :param fixed: list of (fixed) atoms
        :param moving: list of (moving) atoms
        :type fixed,moving: [L{Atom}, L{Atom},...]
        """
        if not len(fixed) == len(moving):
            raise PDBException("Fixed and moving atom lists differ in size")
        length = len(fixed)
        fixed_coord = np.zeros((length, 3))
        moving_coord = np.zeros((length, 3))
        for i in range(0, length):
            fixed_coord[i] = fixed[i].get_coord()
            moving_coord[i] = moving[i].get_coord()
------- sup = SVDSuperimposer()   -------------------------> HERE !!!!!!!!           
------- sup.set(fixed_coord, moving_coord)    --------------> HERE !!!!!!!!
        sup.run()
        self.rms = sup.get_rms()
        self.rotran = sup.get_rotran()

以及https://github.com/biopython/biopython/blob/master/Bio/SVDSuperimposer/init.py#L123

def set(self, reference_coords, coords):
        """Set the coordinates to be superimposed.
        coords will be put on top of reference_coords.
        - reference_coords: an NxDIM array
        - coords: an NxDIM array
        DIM is the dimension of the points, N is the number
        of points to be superimposed.
        """
        # clear everything from previous runs
        self._clear()
        # store coordinates
        self.reference_coords = reference_coords
-----   self.coords = coords    ---------------------------ADD HERE
        n = reference_coords.shape
        m = coords.shape
        if n != m or not (n[1] == m[1] == 3):
            raise Exception("Coordinate number/dimension mismatch.")
        self.n = n[0]

只需在此处添加以下内容：

if array_equal(self.reference_coords, self.coords):
            
            
            raise Exception("Sorry, you are aligning same structure")

需要from numpy import array_equal
不知道SVDSuperimposer（super_imposer.rotran）如何返回
旋转和从坐标的平移馈送到
算法，但可能是你可以让他们甚至从相同的坐标开始
参考和样品？？？我不确定我的只是个假设。

赞(0）回复(0）举报 2023-05-05

我来回答

python 为什么这个源代码给出了不正确的RMSD值？

1条答案

相关问题

热门标签

最新问答