OpenStructure
renumber.py
Go to the documentation of this file.
1 from ost import io, seq, mol, conop
2 from ost import *
3 
4 def Renumber(seq_handle, sequence_number_with_attached_view=1):
5  """
6  Function to renumber an entity according to an alignment between the model sequence
7  and the full-length target sequence. The aligned model sequence or the alignment itself
8  with an attached view needs to be provided. Upon succcess, the renumbered entity is returned.
9 
10  .. code-block:: python
11 
12  from ost.seq.alg import renumber
13  from ost.bindings.clustalw import *
14  ent=io.LoadPDB("path_to_model")
15  s=io.LoadSequence("path_to_full_length_fasta_seqeunce")
16  pdb_seq=seq.SequenceFromChain("model", ent.chains[0])
17  aln=ClustalW(s,pdb_seq)
18  aln.AttachView(1,ent.chains[0].Select(""))
19  e=Renumber(aln.GetSequence(sequence_number_with_attached_view))
20  io.SavePDB(e, "renum.pdb")
21 
22  """
23  if isinstance(seq_handle, seq.SequenceHandle):
24  if seq_handle.HasAttachedView()==False:
25  raise RuntimeError, "Sequence Handle has no attached view"
26  changed_residue_count=0
27  renumberingFlag = False
28  ent_n=mol.CreateEntity()
29  ed=ent_n.EditXCS()
30  c=ed.InsertChain(" ")
31  for pos in range(len(seq_handle)):
32  if seq_handle[pos]!='-':
33  r=seq_handle.GetResidue(pos)
34  if r.IsValid():
35  #print seq_handle[pos],r.number.num,pos+1
36  if r.number.num!=pos+1:
37  changed_residue_count+=1
38  renumberingFlag = True
39  r_n=ed.AppendResidue(c,r.name, mol.ResNum(pos+1))
40  for atom in r.atoms:
41  ed.InsertAtom(r_n,atom.name,atom.pos,atom.prop)
42  else:
43  err='Error: renumbering failed at position %s' %pos
44  raise RuntimeError, err
45  if renumberingFlag == True:
46  err = 'Warning: %s residues have been renumbered!' %changed_residue_count
47  LogInfo(err)
48  conop.ConnectAll(ent_n)
49  return ent_n
50 
51  elif isinstance(seq_handle, seq.AlignmentHandle):
52  if seq_handle.GetSequence(sequence_number_with_attached_view).HasAttachedView()==False:
53  raise RuntimeError, "Sequence Handle has no attached view"
54  dir(seq_handle)
55  counter=0
56  changed_residue_count=0
57  renumberingFlag = False
58  ent_n=mol.CreateEntity()
59  ed=ent_n.EditXCS()
60  c=ed.InsertChain(seq_handle.GetSequence(sequence_number_with_attached_view).GetAttachedView().chains[0].name)
61  for col in seq_handle:
62  if col[0]!='-' and col[1]!='-':
63  if col[0]==col[1]:
64  rnum=seq_handle.GetSequence(sequence_number_with_attached_view).GetResidueIndex(counter)
65  r=seq_handle.GetSequence(sequence_number_with_attached_view).GetResidue(counter)
66  if r.IsValid():
67  if r.number.num!=counter+1:
68  changed_residue_count+=1
69  renumberingFlag = True
70  r_n=ed.AppendResidue(c,r.name, mol.ResNum(counter+1))
71  for atom in r.atoms:
72  ed.InsertAtom(r_n,atom.name,atom.pos,atom.element, atom.b_factor,
73  atom.occupancy, atom.is_hetatom)
74 
75  else:
76  raise RuntimeError("invalide residue at postion %s (renumbering failed)" %(counter))
77  else:
78  raise RuntimeError("residue mismatch at position %d (%s vs %s) (renumbering failed)"%(counter, col[0],col[1]))
79  counter+=1
80  if renumberingFlag == True:
81  err = 'Warning: %s residues have been renumbered!' %changed_residue_count
82  LogInfo(err)
83  conop.ConnectAll(ent_n)
84  return ent_n
85