Thursday, October 30, 2008

interpolate missing residues in pdb file

1: download modeller(http://salilab.org/modeller/wiki/Missing%20residues)
2: you can get FASTA file from RCSB PDB(sequence details)
3: use sequence.py to generate .ali files
#################################
import sys
from modeller import *
from modeller.scripts import complete_pdb

# Get the sequence of the 1qg8 PDB file, and write to an alignment file
chain_id = sys.argv[1]

env = environ()
env.io.atom_files_directory = ['.', '../atom_files']
env.libs.topology.read(file='$(LIB)/top_heav.lib')
env.libs.parameters.read(file='$(LIB)/par.lib')

aln = alignment(env)
mdl = model(env)

# read sequence from pdb list
code = 'XXXX'
mdl.read(file=code, model_segment=('FIRST:' + chain_id, 'LAST:' + chain_id))

# add the sequence to the alignment
aln.append_model(mdl, align_codes=code, atom_files=code)

# read complete pdb with missing residue
fasta='XXXX_' + chain_id + '.fasta.txt'
aln.append(file=fasta, alignment_format='FASTA')

# align them by sequence
aln.malign(gap_penalties_1d=(-500, -300))
aln.write(file='chain'+chain_id+'.ali')
#################################

4: use patch.py to add missing residues.

#################################
# insert the missing residues

import sys
from modeller import *
from modeller.automodel import * # Load the automodel class

log.verbose()
env = environ()

# directories for input atom files
env.io.atom_files_directory = ['.', '../atom_files']

chain_id = sys.argv[1]


class MyModel(automodel):
def select_atoms(self):
if chain_id == 'A' :
return selection(self.residue_range('817', '871'))
elif chain_id == 'C' :
return selection(self.residue_range('1', '7'),
self.residue_range('42', '61'),
self.residue_range('424', '431'))
elif chain_id == 'D' :
return selection(self.residue_range('1', '9'))
elif chain_id == 'E' :
return selection(self.residue_range('1', '8'),
self.residue_range('74', '76'))

a = MyModel(env, alnfile = 'chain' + chain_id + '.ali',
knowns = 'chain' + chain_id, sequence = 'chain'+chain_id+'_fill')
a.starting_model= 1
a.ending_model = 1

a.make()
#################################

No comments: