Hello, all
I have written a script which extracts certain residues from a PDB file, and
writes them out to a new PDB-file.
The script looks like this:
from rdkit import Chem
def create_sub_mol(mol_in, mol_in_conf, bs_atom_list):
"""
Creates a mol of given atoms atoms, e.g. the atoms where deltaSASA is larger
than zero.
"""
empty_mol = Chem.Mol() #create empty mol
mol_rw = Chem.RWMol(empty_mol) #make empty mol editable. This will be our output
conformer = Chem.Conformer(len(bs_atom_list)) # create conformer for new mol.
Otherwise, unable to add xyz
for idx in bs_atom_list:
#get original atom info
original_atom = mol_in.GetAtomWithIdx(idx)
res_info = original_atom.GetPDBResidueInfo()
original_element = original_atom.GetAtomicNum()
original_res_name = res_info.GetResidueName()
original_xyz = list(mol_in_conf.GetAtomPosition(idx)) #
#set new atom, info
new_xyz = Chem.rdGeometry.Point3D(original_xyz[0], original_xyz[1],
original_xyz[2]) #create Point3D object which can be added to conformer
new_atom_idx = mol_rw.AddAtom(Chem.Atom(original_element)) #add atom, returns
new atom idx
print(new_atom_idx)
new_atom = mol_rw.GetAtomWithIdx(new_atom_idx)
conformer.SetAtomPosition(new_atom_idx, new_xyz) #assign new xyz to new mol
new_res_inf = Chem.AtomPDBResidueInfo()
new_res_inf.SetResidueName(original_res_name)
new_atom.SetMonomerInfo(new_res_inf) #invoking this line creates trouble
mol_out = mol_rw.GetMol()
mol_out_conf = mol_out.AddConformer(conformer, assignId = True) #merges
conformer with mol
return mol_out, mol_out_conf
rna = Chem.rdmolfiles.MolFromPDBFile("./1aju.pdb")
rna_conf = rna.GetConformer()
short_list = [65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99,
100, 101, 102, 103, 104, 105, 106, 127, 128, 129, 130, 131, 132, 133, 134, 135,
136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151,
152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167,
168, 169, 170, 171, 323, 324, 325, 326, 327, 328, 329, 330, 331, 332, 333, 334,
335, 336, 337, 338, 339, 340, 341, 342, 343, 344, 345, 346, 347, 348, 349, 350,
351, 352, 353, 354, 355, 356, 357, 358, 359, 360, 361, 362, 363, 364, 365, 366,
367, 368, 369, 370, 371, 372, 373, 374, 375, 376, 377, 378, 379, 380, 381, 382]
rna_extract, rna_extract_conf = create_sub_mol(rna, rna_conf, short_list)
Chem.rdmolfiles.MolToPDBFile(rna_extract, "1aju_extr.pdb")
When I invoke the line with the function SetMonomerInfo() for the new atom, the
coordinates for the extracted residues suddenly get jumbled and compressed, and
make no sense, as opposed to when I create this mol-object without invoking
this line.
What is happening here that I do not understand?
Looking forward to hearing from you all,
Illimar Rekand
Ph.D. candidate,
Brenk-lab, Haug-lab
Department of Biomedicine
Department of Chemistry
University of Bergen
_______________________________________________
Rdkit-discuss mailing list
Rdkit-discuss@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
