--ChemShiftMapFromTwoPeakLists -- script takes two peaklists as input (user edits their ID numbers below) -- for each peak with same peak id in both peaklists it -- calculates the combined chemical shift difference (CCSD) (weighting is defined below) -- writes CCSD of each peak to termina windowl and -- calculates the ave +- stdev of CCSD --for all peaks whose CCSD > STDEV(CCSD) -- it writes a line in MOLMOLs short selection format to select the CA atoms of the residues correspondingto the peaknumbers -- (peaklists should therefore be generated by LUA script PeakListNumberedByResidue) -- it also writes a line for each peak with a given color (edited by user below) -- You can use these selections as follows in Molmol: -- 1) select your molecule by clicking the menu button [Mol] , in Molecule dialog click and click (or click-drag) on desired molecule(s) in list -- 2) Run the user macro ribbons by clicking the menu button [ribbon] to generate a ribbon view of your molecule -- 3) Select menu Primitive-Ribbon-Paint , and inside this menu click button -- 4) alternatively you can type "PaintRibbon atom" on Molmols command line -- 5) copypaste the atom selections output of this Lua script into the Atom field of Molmols selection dialog -- 6) click on [Color] menu to open Molmols Color dialog, select a color and then click on [Atom] --version 1 -F. Damberger 2015-09-18 -- user editable parameters -------------------------------------------------------------------------- ProjectName="myprojectname" PeakListId1 = 1 -- apo PeakListId2 = 2 -- + ligand PeakColor = 2 -- enter an integer corresponding to peakcolors you want to output molmol selection for -- chemical shift scaling Scale1=1 -- weighting of direct dimensions chemical shift (usually =1) Scale2=0.14 -- weighting of indirect dimension chemical shift (typically 0.14 for 15N amide shifts) StdevThreshhold = 1.0 -- 1.0 = molmol output selects residues whose peaks have shifts > 1 x standard deviation -- molmol output parameters MolNum = "1" -- molecule number for Molmol output AtomName = "CA" -- aotmname of atoms to select in Molmol output -- format for output of combined chemical shift changes FormatShift = "%3.3f" -- %N.Mf means format as floating point with M values after the decimal -- end of user editable section ---------------------------------------------------------------------- -- get project and peaklists and check that they are valid choices Project=cara:getProject(ProjectName) if not Project then Project=cara:getProject() end -- try to get project (assuming repository contains only one project) if not Project then error("Please edit top of script to define a ProjectName so it is an existing projectname") end PeakList1=Project:getPeakList( PeakListId1 ) -- apo PeakList2=Project:getPeakList( PeakListId2 ) -- + ligand if not PeakList1 then error("Please edit top of script to define PeakListId1 = Id of an existing peaklist") end if not PeakList2 then error("Please edit top of script to define PeakListId2 = Id of an existing peaklist") end --define functions function square( Val ) return Val*Val end --main script print("\n") print("Output of ChemShiftMapeFromTwoPeakLists\n") print("Selected peaklists:", PeakList1:getId()..":"..PeakList1:getName(), PeakList2:getId()..":"..PeakList2:getName() ) print("Threshold cutoff for chemical shift mapping: "..StdevThreshhold.." times Standard Deviation of combined shift changes" ) print("Color of peaks to select:", PeakColor ) print("\n") ChemShiftMapTable = {} n = 0 -- number of peaks Sum = 0 -- of combined chemical shift changes SqSum = 0-- of combined chemical shift changes Peaks1=PeakList1:getPeaks() CStable = {} -- table of combined chemical shift changes (index is PeakId) PeakColorTable = {} PeakIdTable = {} -- Table of PeakId s (this will get sorted in ascending PeakId order) --determine variance for PeakId1,Peak1 in pairs( Peaks1 ) do Peak2=PeakList2:getPeak(PeakId1) if Peak2 then xy1={} --X=xy1[1],Y=xy1[2] coordinate of Peak1 xy2={} --X=xy2[1],Y=xy2[2] coordinate of Peak2 xy1[1],xy1[2]=Peak1:getPos() xy2[1],xy2[2]=Peak2:getPos() PeakColor2=Peak2:getColor() ShiftDiffSquare= square( Scale1 * ( xy1[1] - xy2[1] ) ) + square( Scale2 * ( xy2[2] - xy2[2] ) ) ShiftDiff = math.sqrt( ShiftDiffSquare ) n = n + 1 PeakIdTable[ n ] = PeakId1 CStable[ PeakId1 ] = ShiftDiff PeakColorTable[ PeakId1 ] = PeakColor2 Sum = Sum + ShiftDiff SqSum = SqSum + ShiftDiffSquare end end Ave=Sum/n Stdev=math.sqrt( ( SqSum - Sum*Sum/n ) / n ) --output result of Standard Deviation calculation --format values Sum=string.format( FormatShift, Sum) SqSum=string.format( FormatShift, SqSum) Ave=string.format( FormatShift, Ave) Stdev=string.format( FormatShift, Stdev ) --print out values print("Number of peaks,Calculated Sum, Square Sum, Average and Variance of combined chemical shifts") print("n", n, "Sum", Sum, "SqSum", SqSum, "Ave", Ave, "Stdev", Stdev ) print("\n") table.sort( PeakIdTable ) -- sort PeakIdTable by PeakId number (before this sort, order is random) print("table of combined chemical shift changes") print("------------------------------------------------------------------------") print("Index","PeakId","Combined ChemShift") for i=1,n do print( i, PeakIdTable[i], string.format(FormatShift, CStable[ PeakIdTable[i] ] ) ) end -- find peaks with cutoff above variance -- and peakcolor = to selected color -- and concatanate corresponding residue numbers for molmol selection PeakListCS = false PeakListColor = false for i=1,n do if CStable[ PeakIdTable[ i ] ] > Stdev*StdevThreshhold then if not PeakListCS then PeakListCS = PeakIdTable[ i ] else PeakListCS = PeakListCS..","..PeakIdTable[ i ] end end if PeakColorTable[ PeakIdTable[ i ] ] == PeakColor then if not PeakListColor then PeakListColor = PeakIdTable[ i ] else PeakListColor = PeakListColor..","..PeakIdTable[ i ] end end end --generate molmol output print("\n") print("Molmol selections: copypaste into Molmol selection window field Atoms") print("Select all ribbon primitives and color by Atom") print("-------------------------------------------------------------------------------") print("Molmol selections for residues with combined shift differences > "..StdevThreshhold.." * stdev" ) print( "#"..MolNum..":"..PeakListCS.."@"..AtomName.."\n") print("Molmol selection for residues whose peaks have peakcolor "..PeakColor) print( "#"..MolNum..":"..PeakListColor.."@"..AtomName.."\n") print( "ChemShiftMapFromTwoPeakLists is done. See terminal window for output." )