-- script to report residues with large deviations from database chemical shifts -- the deviation of assigned chemical shifts from the average values from the database -- are calculated, normalized to the standard deviations from database and -- the square of these values -- averaged over all assigned spins in each residue -- the squareroot of this value "AveDev" is written out for each residue. -- a large value indicates significant deviation from database statistics -- e.g. AveDev > 2.0 should definitely be looked at. -- Large deviations in H atoms may be explained by nearby aromatic rings -- To investigate compare individual spins shifts with database use -- AssignmentReport.lua -- F.Damberger July 1 2004. -- create array for script variables t = {} -- ---------------------------------------------------------------------- --1. Get ProjectName local ProjectNames = {} i = 0 for a,b in pairs(cara:getProjects()) do i = i + 1 ProjectNames[ i ] = b:getName() end t.ProjectName=dlg.getSymbol("Select Project","", unpack( ProjectNames ) ) t.P = cara:getProject( t.ProjectName ) --2. Get output format t.OutputMethod = dlg.getSymbol("Select Output Method","", "output to terminal window", "output to file" ) --3. Get output Filename (if needed) if t.OutputMethod == "output to file" then t.Filename = dlg.getText("Enter the output filename", "", "SideChainDeviations_for_"..t.ProjectName..".txt" ) end --4. Get AtomTypes to Report t.AtomChoice = dlg.getSymbol( "Select which Atoms to check","","H atoms only","ALL" ) -- --------------------------------------------------------------------- -- cycle through all residues in Sequence t.Seq = t.P:getSequence() i = 0 t.Lines = {} for ResNum,Res in pairs( t.Seq ) do t.Assigned = " " t.Unassigned = " " t.Res = t.P:getResidue( ResNum ) t.Typ = t.Res:getType() t.Name = t.Typ:getName() t.Atoms = t.Typ:getAtoms() t.Sys = t.Res:getSystem() t.Sum = 0 t.NumAtoms = 0 -- print("--- "..Res:getType():getName()..ResNum.."---") for AtomId,Atom in pairs( t.Atoms ) do t.ReportThisAtom = false if t.AtomChoice == "H atoms only" then if Atom:getAtomType() == "H" then t.ReportThisAtom = true end elseif t.AtomChoice == "ALL" then if Atom:getAtomType() ~= "O" and Atom:getAtomType() ~= "S" then t.ReportThisAtom = true end else t.ReportThisAtom = false end -- if Atom:getAtomType() ~= "O" and Atom:getAtomType() ~= "S" then if t.ReportThisAtom then t.AtomName = Atom:getName() if Atom:getGroup() ~= nil then -- get GroupName (PseudoAtoms) t.GroupName = Atom:getGroup():getName() else t.GroupName = nil end for SpinId,Spin in pairs( t.Sys:getSpins() ) do if t.AtomName == Spin:getLabel() or t.GroupName == Spin:getLabel() then SpinAssigned = true -- print( "Atom "..Atom:getName().." assigned to "..Spin:getLabel() ) -- compare chemical shifts & write into t.Outliers if needed t.LibraryAve,t.LibrarySd = t.Res:getValue( Atom:getName() ) t.Shift = Spin:getShift() t.Dev = t.Shift - t.LibraryAve t.DevNorm = t.Dev / ( t.LibrarySd / 4 ) t.Sum = t.Sum + ( t.DevNorm * t.DevNorm ) t.NumAtoms = t.NumAtoms + 1 end -- if spin matches atom, then SpinAssigned end -- searched through all spins in system end -- ignore Oxygen & Serine (not assignable) end -- scan through all atoms in residue if t.NumAtoms > 0 then t.Sum = t.Sum / t.NumAtoms t.Sum = math.sqrt( t.Sum ) t.Sum = string.format("%5.3f", t.Sum) -- t.AveDev[ ResNum ] = t.Sum i = i + 1 t.Lines[ i ] = ResNum.." "..t.Sum.." "..t.NumAtoms end end -- scan through entire sequence -- Output the results if t.OutputMethod == "output to terminal window" then for i=1,table.getn( t.Lines ) do print(t.Lines[ i ]) end print( "\n".."Printed out report to terminal window" ) else -- open file to write outfile = io.output( t.Filename ) for i=1,table.getn( t.Lines ) do outfile:write( t.Lines[ i ].."\n" ) end outfile:close() print( "Wrote out report to file "..t.Filename ) end t = nil