#!/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 NormalVector(v): return v/sqrt(v.dotProduct(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[0] ##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): # print filename sse_dict = pickout(sse_annot) # print sse_dict for key in sse_dict.keys(): # print key if key[0] in filename: # print key[0] sse_list=[] chain1 = protein[0] for x in sse_dict[key]: sse_list.append(chain1[x[0]:x[1]]) # print x[0],x[1] # print sse_list return Collection(sse_list) ''' chain1 = protein[0] helix1 = chain1[48:78] DomainA = Collection(chain1[0:36], chain1[123:243]) DomainN = chain1[359:604] DomainP = Collection(chain1[329:359], chain1[604:737]) DomainPRed = chain1[604:737] Loop67 = chain1[811:822] HeliceP= chain1[604:617] helix1 = chain1[48:78] helix2 = chain1[88:121] helix3 = chain1[247:274] helix4 = chain1[288:329] helix4bas = chain1[288:306] helix4haut = chain1[310:329] helix5 = chain1[739:779] helix6 = chain1[788:808] helix7 = chain1[830:852] helix8 = chain1[895:914] helix9 = chain1[930:949] helix10 = chain1[965:984] TM = Collection(helix1,helix2,helix3,helix4,helix5,helix6,helix7,helix8,helix9,helix10)''' #helix=TM #suffix='TM' ## helix is changed to sse_collection def calculatingthings(sse_collection,suffix, motion_type): cm, mi = sse_collection.centerAndMomentOfInertia() ev, axes = mi.diagonalization() print argmin(ev), argmax(ev) axes=axes.inverse() sse_axis = axes[argmin(ev)].asVector().normal() print sse_axis vector = ParticleVector(universe) # print vector graphics = VMD scene = graphics.Scene(scale=1./Units.Ang) if motion_type=="rotation": for atom in sse_collection.atomList(): vector[atom] = sse_axis.cross(atom.position()-cm) print atom.position()-cm print sse_axis.cross(atom.position()-cm) 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)) scene.view() vector=NormalVector(vector) p = [] for i in range(0, len(modes)): mode = NormalVector(modes[i]) projection = mode.dotProduct(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": 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)) scene.view() vector=NormalVector(vector) p = [] pacc = 0 for i in range(6,len(modes)): mode=NormalVector(modes[i]) projection=mode.dotProduct(vector) p.append((i+1, projection**2)) p=array(p) pacc=add.accumulate(p) writeArray(p, 'translation_%s.pics.dat' %suffix) pacc[:,0] = p[:,0] writeArray(pacc, 'translation_%s.esc.dat' %suffix) if motion_type == "indiv_translation": #print "7", len(sse_collection[7]) for counter in xrange(len(sse_collection)): element=sse_collection[counter] print counter, len(sse_collection[counter]) element_cm, element_mi = element.centerAndMomentOfInertia() print element_mi ele_ev, ele_axes = element_mi.diagonalization() ele_axis = ele_axes[argmin(ele_ev)].asVector().normal() indiv_trans_vector = ParticleVector(universe) p=[] pacc=0 for atom in element.atomList(): indiv_trans_vector[atom] = ele_axis indiv_trans_vector=NormalVector(indiv_trans_vector) for i in range(6, len(modes)): mode = NormalVector(modes[i]) projection = mode.dotProduct(indiv_trans_vector) p.append((i+1,projection**2)) p=array(p) writeArray(p, 'indiv_translation_%s_%s.pics.dat'% (suffix, counter+1)) pacc = add.accumulate(p) pacc[:,0] = p[:,0] writeArray(pacc, 'indiv_translation_%s_%s.esc.dat'%(suffix,counter+1)) # p = [] return #print alpha_collection #print alpha_collection.centerAndMomentOfInertia() alpha_collection = make_collection(ALPHA_annot, filename) calculatingthings(alpha_collection, sys.argv[2]+"_ALPHA_CA", sys.argv[3]) beta_collection = make_collection(BETA_annot, filename) calculatingthings(beta_collection, sys.argv[2]+"_BETA_CA", sys.argv[3]) calculatingthings(protein, sys.argv[2]+"_WHOLE_CA",sys.argv[3])