-- CompareProjects.lua -- This script compares the assignments of the corresponding atoms in two projects. -- it reports when the Atom is assigned in only one project or only in the other -- it reports when the Atoms shifts differ by more than a maximum offset (see table below) t = {} -- defines a global table for the script -- USER INSTRUCTIONS: -- 1) Use "Renumber from here" in the Sequence-explorer of the Target project to assign the residues ChainNumbers that match the Residue Number in the SourceProject -- if you have an offset between the sequence numbering in the Source and the ChainNumbers in the Target projects for the corresponding residues, do the following: -- Edit the offset below to add a fixed value onto the Residue number in SourceProject that is used to search for a Residue with the corresponding ChainNr in the TargetProject -- I.e. Residue nr N in SourceProject will be matched to the Residue with ChainNr ( N - NrOffset ) NrOffset = 0 Verbose = false -- set to true to ge an alternative output that is more complete (reports also residues with missing assignments in one project) -- PPM offsets between projects, I did not yet add this to the user options in the dialog, -- The user must edit the lines below to introduce an offset between the chemical shifts -- Chemical shift difference = SourceProject - TargetProject + PpmOffset[ AtomType ] -- i.e. if TargetProject 1H shifts are systematically 0.03 ppm higher than SourceProject 1H shifts, then set PpmOffset["H"] = 0.03 -- Ofcourse, it would be better to calibrate both projects properly before comparison which makes the corrections below unnecessary PpmOffset = {} PpmOffset[ "H" ] = 0.00 -- to correct for offset, set to negative of average difference proj1-proj2 PpmOffset[ "C" ] = 0.00 PpmOffset[ "N" ] = 0.00 -- Do not edit below this line -- Function definitions function CreateProjectNamesComboBox( SelectedProjectName ) ProjectNamesComboBox = gui.createComboBox( t.frm ) SelectedItemIndex = nil Count = 0 for Id,Project in pairs( cara:getProjects() ) do Count = Count + 1 ItemIndex = ProjectNamesComboBox:addItem( Project:getName() ) if Count == 1 then FirstItemIndex = ItemIndex end ProjectNamesComboBox:setCurrentItem( ItemIndex ) if SelectedProjectName == ProjectNamesComboBox:getCurrentText() then SelectedItemIndex = ItemIndex end end -- for all projects if SelectedItemIndex then -- set to previous choice if it exists ProjectNamesComboBox:setCurrentItem( SelectedItemIndex ) return ProjectNamesComboBox else ProjectNamesComboBox:setCurrentItem( FirstItemIndex ) return ProjectNamesComboBox end end function GetSpinWithLabel( Sys, AtName ) if not Sys then return nil end for SpinId,Spin in pairs( Sys:getSpins() ) do if Spin:getLabel() == AtName then return Spin end end -- for all Spins in Sys return nil end -- function GetspinWithLabel function TableSize( Table ) local Size = 0 for IdNum,Entry in pairs( Table ) do Size = Size + 1 end -- for return Size end -- function TableSize function GetResWithNr( Project, Nr ) local Seq = Project:getSequence() if not Seq then error( "Project "..Project:getName().." has no sequence! Aborting CompareProjects.lua ") end for Id,Res in pairs( Seq ) do if Res:getNr() == Nr then return Res end end print("WARNING: function GetResWithNr failed to find Residue in Project "..Project:getName().." with Nr "..Nr ) return nil end -- ----------------------------------------------------------------------------- -- ============Set up menu window for user preferences=================== --0. Check whether there are any projects available ProjectsTable = cara:getProjects() if TableSize( ProjectsTable ) == 0 then error("No projects found in repository") end --4. Create Table of LineEdit Labels for DeltaAtoms t.DeltaAtomLabel = {} t.DeltaAtomLineEdit = {} DeltaAtomTable = {"H","C","N"} DefaultDeltaAtom = {} DefaultDeltaAtom["H"]=0.025 DefaultDeltaAtom["C"]=0.3 DefaultDeltaAtom["N"]=0.3 --1. Create main menu window v = gui.createMainWindow() v:setCaption( "AssignmentReport Preferences" ) t.frm = gui.createGrid( v, 2, false ) v:setCentralWidget( t.frm ) v:show() t.frm:show() -- SourceProject Selector t.SourceProjectListLabel = gui.createLabel( t.frm, "Select Source Project( comparison assignments" ) t.SourceProjectListCB = CreateProjectNamesComboBox( SourceProjectName ) -- TargetProject Selector t.TargetProjectListLabel = gui.createLabel( t.frm, "Select Target Project( reference assignments" ) t.TargetProjectListCB = CreateProjectNamesComboBox( TargetProjectName ) -- Display SourceProjectList Combobox t.SourceProjectListLabel:show() t.SourceProjectListCB:show() -- Display TargetProjectList Combobox t.TargetProjectListLabel:show() t.TargetProjectListCB:show() -- Create ReportOnylResWithDiff Checkbox t.ReportOnlyResWithDiffChBxLabel = gui.createLabel( t.frm, "Report only Residues with differences" ) t.ReportOnlyResWithDiffChBx = gui.createCheckBox( t.frm ) if CompareProjects_ReportOnlyResWithDiff then t.ReportOnlyResWithDiffChBx:setChecked() end t.ReportOnlyResWithDiffChBxLabel:show() t.ReportOnlyResWithDiffChBx:show() -- Create DeltaAtom LineEdits for IdDeltaAtom,DeltaAtom in pairs( DeltaAtomTable ) do t.DeltaAtomLabel[ DeltaAtom ] = gui.createLabel( t.frm, "Enter minimal "..DeltaAtom.." ppm difference to report shifts" ) t.DeltaAtomLineEdit[ DeltaAtom ] = gui.createLineEdit( t.frm ) if not CompareProjects_DeltaAtom then CompareProjects_DeltaAtom = {} end if CompareProjects_DeltaAtom[ DeltaAtom ] then t.DeltaAtomLineEdit[ DeltaAtom ]:setText( CompareProjects_DeltaAtom[ DeltaAtom ] ) else t.DeltaAtomLineEdit[ DeltaAtom ]:setText( DefaultDeltaAtom[ DeltaAtom ] ) end -- Display Label and LineEdit window t.DeltaAtomLabel[ DeltaAtom ]:show() t.DeltaAtomLineEdit[ DeltaAtom ]:show() end -- for all isotopes --9. OK and Cancel Buttons t.okbutton = gui.createPushButton(t.frm, "OK" ) t.cancelbutton = gui.createPushButton( t.frm, "Cancel" ) t.okbutton:show() t.cancelbutton:show() -- ============Callbacks for menu window =================== -- Define Callbacks for the buttons -- cancel button Callback t.cancelbutton:setCallback( gui.event.Clicked, function (self) print( "canceled CompareProjects.lua") v:close() end ) -- end cancel button Callback -- OK button Callback t.okbutton:setCallback( gui.event.Clicked, function (self) -- =============== Determine User Preferences =========================== t.SourceProject = cara:getProject( t.SourceProjectListCB:getCurrentText() ) SourceProjectName = t.SourceProject:getName() t.P1 = t.SourceProject t.TargetProject = cara:getProject( t.TargetProjectListCB:getCurrentText() ) TargetProjectName = t.TargetProject:getName() t.P2 = t.TargetProject --t.ReportAllAssignedRes = ReportAllAssignedResCB:getCurrentValue() if t.ReportOnlyResWithDiffChBx:isChecked() then CompareProjects_ReportOnlyResWithDiff = true else CompareProjects_ReportOnlyResWithDiff = false end -- Create array for DeltaAtoms if it is missing if not CompareProjects_DeltaAtom then CompareProjects_DeltaAtom = {} else end -- get the DeltaAtom values for IdDeltaAtom,DeltaAtom in pairs( DeltaAtomTable ) do CompareProjects_DeltaAtom[ DeltaAtom ] = t.DeltaAtomLineEdit[ DeltaAtom ]:getText() end -- for IdDeltaAtom,DeltaAtom in pairs( DeltaAtomTable ) do -- t.DeltaAtomLabel[ DeltaAtom ] = gui.createLabel( t.frm, "Enter minimal "..DeltaAtom.." ppm difference to report shifts" ) -- t.DeltaAtom[ DeltaAtom ] = gui.createLineEdit( t.frm ) -- end -- ============== Run analysis =========================================== asstext = string.format("%8s"," ass" ) ASStext = string.format("%8s"," ASS" ) NAtext = string.format("%8s"," NA" ) PpmLabelText = string.format("%20s","Source-Target[ppm]") -- Main part of script --t.P1 = cara:getProject( t.projectname1 ) --t.P2 = cara:getProject( t.projectname2 ) t.Seq1 = t.P1:getSequence() -- initialize counters NumResNotAssigned = 0 NumResAssignedOnly1 = 0 NumResAssignedOnly2 = 0 NumResAssignedBoth = 0 NumMissingIn2Not1 = 0 NumMissingIn1Not2 = 0 TotalPpmDiff = {} NumAtomsAssignedInBothProjects = {} NumDiffShifts = {} -- define the table for AtType,Offset in pairs( CompareProjects_DeltaAtom ) do NumDiffShifts[ AtType ] = 0 TotalPpmDiff[ AtType ] = 0 NumAtomsAssignedInBothProjects[ AtType ] = 0 end for ResId1,Res1 in pairs( t.Seq1 ) do ResDoNotDiffer = true TextForRes = "" ResNameText = string.format("%8s",Res1:getType():getLetter()..ResId1) TextForRes = TextForRes.."==========================" CR = "\n" System1 = Res1:getSystem() Res2 = GetResWithNr( t.P2, ResId1+ NrOffset ) if Res2 then System2 = Res2:getSystem() end if System1 then ReportRes1 = ASStext else ReportRes1 = NAtext end if System2 then ReportRes2 = ASStext else ReportRes2 = NAtext end if not System1 and not System2 then NumResNotAssigned = NumResNotAssigned + 1 end if System1 and not System2 then NumResAssignedOnly1 = NumResAssignedOnly1 + 1 ResDoNotDiffer = false end if not System1 and System2 then NumResAssignedOnly2 = NumResAssignedOnly2 + 1 ResDoNotDiffer = false end if System1 and System2 then NumResAssignedBoth = NumResAssignedBoth + 1 end TextForRes = TextForRes..CR..ResNameText..ReportRes1..ReportRes2..PpmLabelText for AtomName,Atom in pairs( Res1:getType():getAtoms() ) do SpinAssignedInOnlyOne = false ShiftsDiffer = false AtomNameText = string.format("%8s", AtomName ) Spin1 = GetSpinWithLabel( System1, AtomName ) Spin2 = GetSpinWithLabel( System2, AtomName ) if Spin1 then ReportSpin1 = asstext else ReportSpin1 = NAtext end if Spin2 then ReportSpin2 = asstext else ReportSpin2 = NAtext end if Spin1 and Spin2 then -- check if shifts differ significantly PpmDiff = Spin1:getShift() - Spin2:getShift() TotalPpmDiff[ Spin1:getAtomType() ] = TotalPpmDiff[ Spin1:getAtomType() ] + PpmDiff NumAtomsAssignedInBothProjects[ Spin1:getAtomType() ] = NumAtomsAssignedInBothProjects[ Spin1:getAtomType() ] + 1 if math.abs( PpmDiff + PpmOffset[ Spin1:getAtomType() ] ) > tonumber( CompareProjects_DeltaAtom[ Spin1:getAtomType() ] ) then -- prepare report on different shifts Shift1 = string.format( "%8.3f", Spin1:getShift() ) Shift2 = string.format( "%8.3f", Spin2:getShift() ) NumDiffShifts[ Spin1:getAtomType() ] = NumDiffShifts[ Spin1:getAtomType() ] + 1 --Flagged = true ShiftsDiffer = true ShiftSpin1 = Shift1 ShiftSpin2 = Shift2 -- next line testing only if not Verbose then print( ResNameText.." "..AtomNameText..ShiftSpin1..ShiftSpin2..string.format( "%20.3f", PpmDiff ) ) end end end if Spin1 and not Spin2 then --ReportSpin2 = ReportSpin2.."NA" -- in TargetProject only!" SpinAssignedInOnlyOne = true NumMissingIn2Not1 = NumMissingIn2Not1 + 1 end if not Spin1 and Spin2 then --ReportSpin2 = ReportSpin2.."NA" --missing in SourceProject only!" SpinAssignedInOnlyOne = true NumMissingIn1Not2 = NumMissingIn1Not2 + 1 end if System1 or System2 then -- only break down by atom if at least one of the two systems is assigned --if SpinAssignedInOnlyOne then --print this atoms results --print("........ "..AtomName.." .........") if ShiftsDiffer then TextForRes = TextForRes..CR..AtomNameText..ShiftSpin1..ShiftSpin2..string.format( "%20.3f", PpmDiff ) else TextForRes = TextForRes..CR..AtomNameText..ReportSpin1..ReportSpin2 end --end end if SpinAssignedInOnlyOne or ShiftsDiffer then ResDoNotDiffer = false end end if CompareProjects_ReportOnlyResWithDiff == false then if Verbose then print( TextForRes ) end else if ResDoNotDiffer then if Verbose then print( TextForRes ) end end end end function For8( text ) return string.format( "%8s", text ) end function For20( text ) return string.format( "%20s", text ) end print("-----------------------------------------------------") print( For8("Res").." "..For8("Atom")..For8("proj1")..For8("proj2")..For20("diff proj1-proj2") ) print(" ================================") if Verbose then print(" KEY ") print("NA = not assigned, ASS = residue assigned, ass = atom assigned") end print("Shifts of atoms are reported if atom shifts differ by more than:") print(" ================================") for IdDeltaAtom,DeltaAtom in pairs( CompareProjects_DeltaAtom ) do print(IdDeltaAtom..": "..DeltaAtom) end if Verbose then print("If no atom shifts differ then only the Residue line ASS ASS is printed") print("Then ass is replaced by the shift values") end print(" ================================") print(" SUMMARY ") print(" ================================") print("----Assigned residues ----") print( string.format( "%4d", NumResNotAssigned ).." : unassigned in both projects " ) print( string.format( "%4d", NumResAssignedOnly1 ).." : assigned only in Project "..SourceProjectName.." : " ) print( string.format( "%4d", NumResAssignedOnly2 ).." : assigned only in Project "..TargetProjectName.." : " ) print( string.format( "%4d", NumResAssignedBoth ).." : assigned in both projects " ) print("------Total number of atom assignments missing -------") print( string.format( "%4d", NumMissingIn1Not2 ).." : only in Project "..SourceProjectName ) print( string.format( "%4d", NumMissingIn2Not1 ).." : only in Project "..TargetProjectName ) print("---- Total number of atoms assigned in both projects and their average ppm difference SourceProject-TargetProject ---") for AtType,Offset in pairs( CompareProjects_DeltaAtom ) do AvePpmDiff = string.format( "%3.2f", TotalPpmDiff[ AtType ]/NumAtomsAssignedInBothProjects[ AtType ] ) print( "Average ppm difference proj1-proj2 for atomType "..AtType.." : "..AvePpmDiff ) end print("---- Total number of atoms whose shifts differ---") for AtType,Offset in pairs( CompareProjects_DeltaAtom ) do print( string.format( "%4d", NumDiffShifts[ AtType ] ).." : "..AtType.." shifts differing by more than "..Offset ) end print("Finished with comparison of projects: SourceProject "..SourceProjectName.." and TargetProject "..TargetProjectName ) v:close() end ) -- end OK button Callback