#!/opt/local/bin/python2.6 '''Adapted from helix bundle rotation code written by Nathalie for TMM@''' import MMTK from MMTK import * from MMTK.PDB import PDBConfiguration from MMTK.Proteins import Protein from MMTK.NormalModes import NormalModes from numpy import sqrt, add, array, argmax, argsort, argmin from Scientific import * from Scientific.IO.ArrayIO import writeArray from Scientific.Visualization import VMD import sys from rpy import * import PlotTableR, util def massWeightedNormalVector(v): return v/sqrt(v.massWeightedDotProduct(v)) filename=sys.argv[1].replace(".modes","") modes = MMTK.load(sys.argv[1]) #should have the input modes file for the TIM barrel structure in question universe = modes.universe protein = universe.protein protein.normalizeConfiguration('Ir') protein.writeToFile(sys.argv[2]+"_"+sys.argv[3]+"_from_modes.pdb")#sys.argv[2] should be filename, sys.argv[3] should be motion_type ##sys.argv[4] should be mode type ##the structure annotation part ALPHA_annot = open("/Users/sandhyatiwari/Documents/TIM Barrel PDBs/ALPHA_from_pdb.txt","r") BETA_annot = open("/Users/sandhyatiwari/Documents/TIM Barrel PDBs/BETA_from_pdb.txt","r") def pickout(file1): name_list={} name='' first=0 for line in file1: line=line.rstrip() row=line.split("\t") if len(row[0]) ==4: name=row[0] first = int(row[1]) name_list[(name,first)]=[] old_line=[] if len(row[0])<4 and "#" not in line: start=int(row[0])-first end=int(row[1])-first name_list[(name,first)]=old_line old_line.append((start,end)) return name_list def make_collection(sse_annot,filename): sse_dict = pickout(sse_annot) for key in sse_dict.keys(): if key[0] in filename: sse_list=[] chain1 = protein[0] for x in sse_dict[key]: sse_list.append(chain1[x[0]:x[1]]) return Collection(sse_list) ## helix is changed to sse_collection def calculatingthings(sse_collection,suffix, motion_type): cm, mi = sse_collection.centerAndMomentOfInertia() ev, axes = mi.diagonalization() sse_axis = axes[argmax(ev)].asVector().normal() vector = ParticleVector(universe) if motion_type=="rotation": for atom in sse_collection.atomList(): vector[atom] = sse_axis.cross(atom.position()-cm) vector=massWeightedNormalVector(vector) p = [] for i in range(6, len(modes)): mode = massWeightedNormalVector(modes[i]) projection = mode.massWeightedDotProduct(vector) p.append((i+1,projection**2)) p = array(p) writeArray(p, 'rotation_%s.pics.dat' %suffix) pacc = add.accumulate(p) pacc[:,0] = p[:,0] writeArray(pacc,'rotation_%s.esc.dat' %suffix) if motion_type == "translation": graphics = VMD scene = graphics.Scene(scale=1./Units.Ang) for atom in sse_collection.atomList(): vector[atom] = sse_axis scene.addObject(graphics.Arrow(atom.position(),atom.position()+1.*vector[atom],0.01)) scene.addObject(graphics.Cylinder(cm-2.*Units.nm*sse_axis, cm+2.*Units.nm*sse_axis, 0.05)) vector=massWeightedNormalVector(vector) p = [] pacc = 0 for i in range(6,len(modes)): ##changed from 0 mode=massWeightedNormalVector(modes[i]) projection=mode.massWeightedDotProduct(vector) p.append((i+1, projection**2)) p=array(p) pacc=add.accumulate(p) writeArray(p, 'translation_nontriv_raw_%s.pics.dat' %suffix) pacc[:,0] = p[:,0] writeArray(pacc, 'translation_nontriv_raw_%s.esc.dat' %suffix) scene.writeToFile("translation_nontriv_raw_scene_%s.vmd"%suffix) if motion_type == "indiv_translation": graphics = VMD scene = graphics.Scene(scale=1./Units.Ang) for counter in xrange(len(sse_collection)): element=sse_collection[counter] element_cm, element_mi = element.centerAndMomentOfInertia() ele_ev, ele_axes = element_mi.diagonalization() ele_axis = ele_axes[argmin(ele_ev)].asVector().normal() indiv_trans_vector = ParticleVector(universe) for atom in element.atomList(): indiv_trans_vector[atom] = ele_axis scene.addObject(graphics.Arrow(atom.position(),atom.position()+1.*indiv_trans_vector[atom],0.01)) scene.addObject(graphics.Cylinder(cm-2.*Units.nm*sse_axis, cm+2.*Units.nm*sse_axis, 0.05)) indiv_trans_vector=massWeightedNormalVector(indiv_trans_vector) for i in range(6, len(modes)): ## changed from 0 to 6 to ignore trivial mode = massWeightedNormalVector(modes[i]) projection = mode.massWeightedDotProduct(indiv_trans_vector) p.append((i+1,projection**2)) p=array(p) writeArray(p, 'indiv_translation_nontriv_%s_%s.pics.dat'% (suffix, counter+1)) pacc = add.accumulate(p) pacc[:,0] = p[:,0] writeArray(pacc, 'indiv_translation_nontriv_%s_%s.esc.dat'%(suffix,counter+1)) p = [] scene.writeToFile("indiv_translation_nontriv_scene_%s.vmd"% suffix) if motion_type=="indiv_slide": p = [] pacc = 0 graphics = VMD scene = graphics.Scene(scale=1./Units.Ang) for counter in xrange(len(sse_collection)): element=sse_collection[counter] element_cm, element_mi = element.centerAndMomentOfInertia() ele_ev, ele_axes = element_mi.diagonalization() ele_axis = ele_axes[argmax(ele_ev)].asVector().normal() indiv_slide_vector = ParticleVector(universe) indiv_rot_vector = sse_axis.cross(cm - element_cm) indiv_slide_v = sse_axis.cross(indiv_rot_vector) if (element_cm - cm) * indiv_slide_v <0: indiv_slide_v = indiv_slide_v * -1 #indiv_slide_v = massWeightedNormalVector(indiv_slide_v) for atom in element.atomList(): indiv_slide_vector[atom] = indiv_slide_v scene.addObject(graphics.Arrow(atom.position(),atom.position()+1.*indiv_slide_vector[atom],0.01)) scene.addObject(graphics.Cylinder(cm-2.*Units.nm*sse_axis, cm+2.*Units.nm*sse_axis, 0.05)) indiv_slide_vector =massWeightedNormalVector(indiv_slide_vector) for i in range(6,len(modes)): ## changed from 0 to 6 to exclude first mode = massWeightedNormalVector(modes[i]) projection = mode.massWeightedDotProduct(indiv_slide_vector) p.append((i+1,projection**2)) p=array(p) writeArray(p, 'indiv_slide_nontriv_%s_%s.pics.dat'% (suffix, counter+1)) pacc = add.accumulate(p) pacc[:,0] = p[:,0] writeArray(pacc, 'indiv_slide_nontriv_%s_%s.esc.dat'% (suffix, counter+1)) p=[] scene.writeToFile("indiv_slide_scene_nontriv_%s.vmd"% suffix) return if __name__ == '__main__': if "indiv" in sys.argv[3]: alpha_collection = make_collection(ALPHA_annot, filename) calculatingthings(alpha_collection, sys.argv[2]+"_ALPHA_"+sys.argv[4], sys.argv[3]) beta_collection = make_collection(BETA_annot, filename) calculatingthings(beta_collection, sys.argv[2]+"_BETA_"+sys.argv[4], sys.argv[3]) else: alpha_collection = make_collection(ALPHA_annot, filename) calculatingthings(alpha_collection, sys.argv[2]+"_ALPHA_"+sys.argv[4], sys.argv[3]) beta_collection = make_collection(BETA_annot, filename) calculatingthings(beta_collection, sys.argv[2]+"_BETA_"+sys.argv[4], sys.argv[3]) calculatingthings(protein, sys.argv[2]+"_WHOLE_"+sys.argv[4],sys.argv[3])