Python pybel.readfile方法代码示例

# 需要导入模块: import pybel [as 别名]
# 或者: from pybel import readfile [as 别名]
def readfile(format, filename, opt=None, lazy=False):
    if format == 'mol2':
        if __version__ < '2.4.0':
            warnings.warn('OpenBabel 2.3.2 does not support writing data in '
                          'comments ["-xc"]. Please upgrade to OB 2.4')
        if opt:
            opt['c'] = None
        else:
            opt = {'c': None}
    if lazy and format == 'mol2':
        return _filereader_mol2(filename, opt=opt)
    elif lazy and format == 'sdf':
        return _filereader_sdf(filename, opt=opt)
    elif lazy and format == 'pdb':
        return _filereader_pdb(filename, opt=opt)
    else:
        return pybel.readfile(format, filename, opt=opt) 

# 需要导入模块: import pybel [as 别名]
# 或者: from pybel import readfile [as 别名]
def mol_from_file(file_path: str, file_format: str = 'xyz'):
    """
    Args:
        file_path(str)
        file_format(str): allow formats that open babel supports
    """
    mol = [r for r in pybel.readfile(format=file_format,
                                     filename=file_path)][0]
    return mol 

# 需要导入模块: import pybel [as 别名]
# 或者: from pybel import readfile [as 别名]
def _load_pybel(self, input, input_type, **kwargs):
        """
        The internal function to load a molecule using rdkit engine.

        """
        # available variables
        pybel_mol = None
        creator = None

        if input_type == 'xyz':
            if os.path.isfile(input):
                creator = ('XYZ', input)
                gen = pybel.readfile("xyz", input)
                mols = list(gen)
                if len(mols) == 1:
                    pybel_mol = mols[0]
                else:
                    warnings.warn('More than one Molecule object is created and is returned as a generator.')
                    return self._multiple_molecules(mols, creator)
            else:
                msg = "The input '%s' is not a valid XYZ input file."%input
                raise ValueError(msg)

        if pybel_mol is None:
            msg = 'The input is not a legit %s string. Refere to the error message that was already displayed by Pybel library!' % \
              creator[0]
            raise ValueError(msg)

        # if the molecule is already being created warn the user
        if self.pybel_molecule and self.creator:
            msg = "The molecule was already built using %s, is now overwritten by %s" % (
            str(self.creator), str(creator))
            warnings.warn(msg)

        self.creator = creator
        self.pybel_molecule = pybel_mol

        self._init_attributes()

        if self.creator[0]=='XYZ':
            self.to_xyz() 

# 需要导入模块: import pybel [as 别名]
# 或者: from pybel import readfile [as 别名]
def get_ring_from_aa(residueid):

    class AAselect(Select):

        def accept_residue(self, residue):
            # print residue.get_full_id()[3][1],residueid
            if str(residue.get_full_id()[3][1]) == residueid:
                return 1
            else:
                return 0
    ptemp = PDBParser(QUIET=True)  # disable warnings
    stemp = ptemp.get_structure(
        pdbname, projectdir + 'pdbs/' + pdbname + '.pdb')
    temp_aa_id = residueid

    io = PDBIO()
    io.set_structure(stemp)
    io.save(projectdir + 'temp/' + residueid + '.pdb', AAselect())

    mol = pybel.readfile("pdb", projectdir + 'temp/' +
                         residueid + '.pdb').next()
    # print hetflag
    rings = getattr(mol, "OBMol").GetSSSR()
    ringlist = []
    for ring in rings:
        center = Vector(0.0, 0.0, 0.0)
        members = ring.Size()
        if ring.IsAromatic():
            atomlist = []
            atomnames = []
            atomvectors = []
            for atom in mol:
                if ring.IsMember(atom.OBAtom):
                    a_vector = Vector(getattr(atom, 'coords'))
                    center += a_vector
                    atomlist.append(atom.idx)
                    atomvectors.append(a_vector)
                    atomnames.append(getattr(atom, 'type'))
            center = center / members
            normal = center - a_vector  # vector in plane
            normal1 = center - atomvectors[0]
            normal2 = center - atomvectors[2]
            normal = Vector(np.cross([normal1[0],normal1[1],normal1[2]],[normal2[0],normal2[1],normal2[2]]))
            ringlist.append([atomlist, center, normal, atomnames, atomvectors])
    return ringlist 

# 需要导入模块: import pybel [as 别名]
# 或者: from pybel import readfile [as 别名]
def get_hydrogen_from_aa(residueid):

    class AAselect(Select):

        def accept_residue(self, residue):
            # print residue.get_full_id()[3][1],residueid
            if str(residue.get_full_id()[3][1]) == residueid:
                return 1
            else:
                return 0
    ptemp = PDBParser(QUIET=True)
    stemp = ptemp.get_structure(
        pdbname, projectdir + 'pdbs/' + pdbname + '.pdb')
    temp_aa_id = residueid

    io = PDBIO()
    io.set_structure(stemp)
    io.save(projectdir + 'temp/' + residueid + '.pdb', AAselect())

    mol = pybel.readfile("pdb", projectdir + 'temp/' +
                         residueid + '.pdb').next()

    mol.OBMol.AddHydrogens(False, True, 7.4)
    # print hetflag
    donors = []
    for atom in mol:
        if getattr(atom, 'OBAtom').IsHbondDonor():
            chargevector = Vector(getattr(atom, 'coords'))
            # print getattr(atom,'type')," is Donor",chargevector
            temphatoms = []
            for neighbor in pybel.ob.OBAtomAtomIter(atom.OBAtom):
                neighbor = pybel.Atom(neighbor)
                if getattr(neighbor, 'type') == "H":
                    # print "neighbor
                    # Atom",getattr(neighbor,'type'),"Coords:",getattr(neighbor,'coords')
                    temphatoms.append(Vector(getattr(neighbor, 'coords')))

            donors.append([getattr(atom, 'type'), chargevector, temphatoms,getattr(atom, 'OBAtom').IsHbondAcceptor()])

        if getattr(atom, 'OBAtom').IsHbondAcceptor():
            chargevector = Vector(getattr(atom, 'coords'))
            #print getattr(atom, 'type'),chargevector,'acceptor!'

    return donors 

# 需要导入模块: import pybel [as 别名]
# 或者: from pybel import readfile [as 别名]
def voxelization(pdb_path, aux_input_path, r, N):
        # Read the pdb file and the auxilary input file
        ppdb = PandasPdb().read_pdb(pdb_path)
        protein_df = ppdb.df['ATOM']
        content = []
        with open(aux_input_path) as in_strm:
            for line in in_strm.readlines():
                l = line.replace('\n','')
                idx = l.index(':')
                content.append(l[idx+1:None])
        in_strm.close()
        resi =  [int(x) for x in content[0].split(' ')]
        if len(content[1]) != 0:
            pocket_df = protein_df[protein_df['residue_number'].isin(resi)]
            pocket_coords = np.array([pocket_df['x_coord'], pocket_df['y_coord'], pocket_df['z_coord']]).T
            pocket_center = np.array([int(x) for x in content[1].split(' ')])
        else:
            print('No center is provided')
            pocket_df = protein_df[protein_df['residue_number'].isin(resi)]
            pocket_coords = np.array([pocket_df['x_coord'], pocket_df['y_coord'], pocket_df['z_coord']]).T
            pocket_center = np.mean(pocket_coords, axis = 0)
        protein_coords = np.array([protein_df['x_coord'], protein_df['y_coord'], protein_df['z_coord']]).T
        pocket_coords = pocket_coords - pocket_center # center the pocket to 0,0,0
        protein_coords = protein_coords - pocket_center # center the protein according to the pocket center
        inertia = np.cov(pocket_coords.T)
        e_values, e_vectors = np.linalg.eig(inertia)
        sorted_index = np.argsort(e_values)[::-1]
        sorted_vectors = e_vectors[:,sorted_index]
        # Align the first principal axes to the X-axes
        rx = vrrotvec(np.array([1,0,0]),sorted_vectors[:,0])
        mx = vrrotvec2mat(rx)
        pa1 = np.matmul(mx.T,sorted_vectors)
        # Align the second principal axes to the Y-axes
        ry = vrrotvec(np.array([0,1,0]),pa1[:,1])
        my = vrrotvec2mat(ry)
        transformation_matrix = np.matmul(my.T,mx.T)
        # transform the protein coordinates to the center of the pocket and align with the principal
        # axes with the pocket
        transformed_coords = (np.matmul(transformation_matrix,protein_coords.T)).T
        # Generate a new pdb file with transformed coordinates
        ppdb.df['ATOM']['x_coord'] = transformed_coords[:,0]
        ppdb.df['ATOM']['y_coord'] = transformed_coords[:,1]
        ppdb.df['ATOM']['z_coord'] = transformed_coords[:,2]
        output_trans_pdb_path = aux_input_path[0:-4] + '_trans.pdb'
        print('Output the binding pocket aligned pdb file to: ' + output_trans_pdb_path)
        ppdb.to_pdb(output_trans_pdb_path)
        output_trans_mol2_path = output_trans_pdb_path[0:-4] + '.mol2'
        print('Output the binding pocket aligned mol2 file to: ' + output_trans_mol2_path)
        mol = pybel.readfile('pdb',output_trans_pdb_path).next()
        mol.write('mol2',output_trans_mol2_path, overwrite = True)    
        # Grid generation and DFIRE potential calculation
        print('...Generating pocket grid representation')
        pocket_grid = Grid3DBuilder.build(protein_coords, output_trans_mol2_path, r, N)
        print('...Generating pocket voxel representation')
        pocket_voxel = site_voxelization(pocket_grid, N + 1)
        pocket_voxel = np.expand_dims(pocket_voxel, axis = 0)
#        np.save('voxel_rep', pocket_voxel)
        return pocket_voxel