from xtra_wrapper import * import util # # sending of parameters to and from zope in the below code is somewhat ugly # # Be careful when relying on the order of stuff in tables: # it seems that tables get reversed when they are converted to strings in zope/biaz # if reversing those back, be sure to reverse every table that are expected to be in sync as well. # Avoid this if possible. I've got the feeling that we are relying on undefined behaviour with this. # # Also do not rely on the order of stuff in dictionaries, # note that the parameters items and results are dictionaries. # def get_kinkPoint_suggestions(items, results): """ parses arguments from tmm-web interface and makes suggestion for kinkpoint for each helix in the TM bundle """ results=items # #use function from xtrawrapper to retrieve some data #place before conversion of strings to lists # helices, __, residues, names, __ = retrieve_data(items) #convert strings to lists #reversing lists as they get reversed by biaz/zope for k,v in results.items(): #results[k] is a string representing a list if len(v) > 0 and v[0]=="[" and v[len(v)-1] == "]": newlist=util.listString2StringList(v) newlist.reverse() results[k]=newlist #extract protein (this could be done from pdb-file) modefilename=results['modefile_name'] modes=load(modefilename) protein=modes.universe[0] #get suggestsions for the kinkpoint print "searching for kinkpoints" kinkpoints=findKink( helices, residues, protein ) print results['jobid'] # #update parameters to reflect users input # new_tm_h=new_tm_helices(items) #remove whitespace names from names print names for n in range(len(names)): if len(names[n].strip()) == 0: names[n]='' new_h_start=new_helix_start(items) new_h_end=new_helix_end(items) # # Reverse the lists # # biaz seems to give lists to zope in the reverse order of how they got in. # reverse all lists here to get the intended sequences # it might be safer to move all naming responsibilities to here and either discard the helix id entry or forward it in a table as well: # names.reverse() results['h_name'] = names new_h_start.reverse() results['h_start']=util.intList2StringList ( new_h_start ) new_h_end.reverse() results['h_end']=util.intList2StringList ( new_h_end ) kinkpoints.reverse() results['kinkpoints'] = util.intList2StringList( kinkpoints ) new_tm_h.reverse() results['tm_helices'] = util.intList2StringList( new_tm_h ) #for diagnostics ## print "results" ## for k, v in results.items(): ## print (k, v) return results # # # Function for updating paramters as specified by user # # def new_helix_start(items): """ updates the h_start list in accordance with helices """ new_h_start=[0]*len(items['h_start']) for k, v in items.items(): if k[0:12] == 'helix_start_': new_h_start[int(k[12:len(k)])]=v return new_h_start def new_helix_end(items): """ updates the h_start list in accordance with helices """ new_h_end=[0]*len(items['h_end']) for k, v in items.items(): if k[0:10] == 'helix_end_': new_h_end[int(k[10:len(k)])]=v return new_h_end def new_tm_helices(items): """ changes tm_helices to reflect which variables tm_helix_## exists These exist if the user have selected them as tm-helices in the web interface """ new_tm_selection=[0]*len(items['tm_helices']) for k,v in items.items(): if k[0:9] == "tm_helix_": new_tm_selection[int(k[9:len(k)])]=1 return new_tm_selection # # # # # # Functions for suggesting kinkpoint. # If these get computationally expensive, make sure to include the bundle paramter and only compute for those helices in bundle # # The list returned by theese functions should have the same size as the list helices and indicate a kinkpoint at the postitons corresponding to a helix in the bundle # The other positions can be undefined # the kinkpoint is the pdb-residuenumber at the center of the region of kink # # Note: When reading the list helices: first postion of each member in helices is the chain identifier. The helix starts at the second postion def findKink(pdbIndexHelices, residues, protein): """ This function relies heavily on the order in which stuff is constructed. pdbIndexHelices is a list of lists representing the helices in the protein. The list representing a helix has the chain identifier prepended to it residues is a list of the pdbIndecies (residuenumber) to every residue in the protein protein is a MMTK object created from the same pdb-file as is used to construct the system that pdbIndexhelices was extracted from returns a list of pdbIndecies representing kinkpoints for the different helices """ #find make a list of all chain identifiers in the order they appear in pdb-file chains=[] lastChain=None for h in pdbIndexHelices: if h[0] != lastChain: chains.append(h[0]) lastChain=h[0] helices=[] #make MMTK-objects that represent helices for h in pdbIndexHelices: chain_nr = util.get_chain_nr(h[0], chains) chain = protein[chain_nr] start_res = h[1] start_res = util.findIndex(start_res, residues, chain_nr) h_len = len(h)-1 end_res = h[h_len] end_res = util.findIndex(end_res, residues, chain_nr) helices.append(chain[start_res:end_res]) # #find the index in the helix that we believe might be the kinkpoint # kinksInHelix=findKinksInMMTKObject(helices) #convert indecies in helices to residuenumbers kinksResidueNumber=[] for i in range(len(kinksInHelix)): kinksResidueNumber.append(pdbIndexHelices[i][kinksInHelix[i]+1]) # add one because elements of pdbIndexHelices starts with chain identifier return kinksResidueNumber # # # Functions that find kinkpoint in a helix # one of these should be plugged into findKink # # # def kinkAtCenter(helices): kinks=[] for h in helices: kinks.append(h[1] + midPoint(h[1:len(h)])) return kinks def findKinksInMMTKObject(helices): """ returns a list of the 0-based indecies to the residues in the helices at which they are suggested to kink """ kinkpoints=[] for helix in helices: kinkpoints.append( naiveKinkFinder(helix, trim=2) ) #kinkpoints.append( propensityKink(helix, trim=2) ) return kinkpoints # #Aux function for kinkAtCenter # def midPoint(list): """ returns the index to the middle member of the list if the list has even length the closest index below will be used """ middle=len(list)/2 if len(list)%2 ==0: return middle else: return middle+1 # # #Aux functions for findKinksInMMTKObject # #One of these should be plugged into findKinksInMMTK object # # def naiveKinkFinder(helix, trim=0): """ Kinks at Proline closest to center If no proline present: kinks at glycine closest to center If no glycine present: kinks at center If len(helix) < 1+2*trim: kinks at center trim: Ignore first and last n residues, where n=trim """ return propensityKink(helix, trim, windowsize=1, propensities={"Pro":2, "Gly":1} ) # The table kinkpropensities is tweaked to make reasonable suggestions for kinkpoints based on the helix primary structure. # The table is built around data found in Kumar and Bansal 1998 (stastical analysis of structures of globular proteins) # The values are modified a bit to favourise kinks at Proline and Glycine in particular, but also to give residues with low alfa-helix propensities or beta branched sidechains a slight preference for kink # Negative values are assigned to some residues that seems to stabilize the helix against kinking (see Kumar and Bansal) # oldkinkpropensities={"Equ":25.0, "Pro":57.8, "Ala":20.6, "Glu":18.9, "Ile":33.4, "Tyr":36.1, "Trp":35, "Gln":28.1, "His": 18.8, "Gly":46.2} kinkpropensities={"Pro":32.8, "Ala":-4.4, "Glu":-6.1, "Ile":8.4, "Tyr":11.1, "Trp":10, "Gln":3.1, "His": -6.2, "Gly":21.2} # # This one could probably work quit OK. # Need some testing and kinkpropensities need further tweaking # def propensityKink(helix, trim=0, windowsize=4, propensities=kinkpropensities): """ Slides a window of size windowsize along helix and calculates that windows overall kinkpropensity. The middle residue (or the last residue in the lower half) for the highest scoring window is returned Ties are resolved by minimizing distance to center returns the 0-based index of the residue at the kinkpoint trim: the function will ignore first and last trim residues """ #return midpoint for small helices #TM-helices is usually longer than 12 residues if len(helix) <= windowsize + 2*trim: return len(helix)/2 #sum all windows and store the tuple (sum, indexOfKinkPoint) in sums sums=[] for windowstart in range(trim, len(helix)-windowsize-trim): sum=0 for r in helix[windowstart:windowstart+windowsize].sequence(): if r in kinkpropensities.keys(): sum+=kinkpropensities[r] #centroidity = length of helix - residue distance from center centroidity=len(helix) - math.fabs(len(helix)-(windowsize/2)) sums.append( (sum, centroidity, windowstart+windowsize/2) ) sums.sort() __, __, bestIndex = sums[len(sums)-1] return bestIndex