diff --git a/actions/ost-compare-structures b/actions/ost-compare-structures index 7104d26d7be36c9c545da89bf80a1bd056183112..6121d05103fce84047d159a4d485b8a07d1e5a57 100644 --- a/actions/ost-compare-structures +++ b/actions/ost-compare-structures @@ -15,7 +15,7 @@ from ost.io import (LoadPDB, LoadMMCIF, MMCifInfoBioUnit, MMCifInfo, MMCifInfoTransOp, ReadStereoChemicalPropsFile) from ost import PushVerbosityLevel from ost.mol.alg import (qsscoring, Molck, MolckSettings, lDDTSettings, - lDDTScorer) + lDDTScorer, CheckStructure) from ost.conop import CompoundLib from ost.seq.alg.renumber import Renumber @@ -407,14 +407,34 @@ def _Main(): # Read the input files ost.LogInfo("Reading model from %s" % opts.model) models = _ReadStructureFile(opts.model) - if opts.molck: - for model in models: - _MolckEntity(model, opts) ost.LogInfo("Reading reference from %s" % opts.reference) references = _ReadStructureFile(opts.reference) if opts.molck: + for i in range(len(references)): + _MolckEntity(references[i], opts) + references[i] = references[i].CreateFullView() + for i in range(len(models)): + _MolckEntity(models[i], opts) + models[i] = models[i].CreateFullView() + if opts.structural_checks: + stereochemical_parameters = ReadStereoChemicalPropsFile( + opts.parameter_file) + ost.LogInfo("Performing structural checks for reference(s)") for reference in references: - _MolckEntity(reference, opts) + CheckStructure(reference, + stereochemical_parameters.bond_table, + stereochemical_parameters.angle_table, + stereochemical_parameters.nonbonded_table, + opts.bond_tolerance, + opts.angle_tolerance) + ost.LogInfo("Performing structural checks for model(s)") + for model in models: + CheckStructure(model, + stereochemical_parameters.bond_table, + stereochemical_parameters.angle_table, + stereochemical_parameters.nonbonded_table, + opts.bond_tolerance, + opts.angle_tolerance) if len(models) > 1 or len(references) > 1: ost.LogInfo( @@ -479,11 +499,9 @@ def _Main(): radius=opts.inclusion_radius, sequence_separation=opts.sequence_separation, sel=opts.selection, - structural_checks=opts.structural_checks, + structural_checks=False, consistency_checks=opts.consistency_checks, label="lddt") - stereochemical_parameters = ReadStereoChemicalPropsFile( - opts.parameter_file) if opts.verbosity > 3: lddt_settings.PrintParameters() # Perform single chain scoring @@ -492,16 +510,15 @@ def _Main(): for aln in qs_scorer.alignments: # Get chains and renumber according to alignment (for lDDT) ch_ref = aln.GetSequence(0).GetName() - reference = Renumber(aln.GetSequence(0)) + reference = Renumber(aln.GetSequence(0)).CreateFullView() ch_mdl = aln.GetSequence(1).GetName() - model = Renumber(aln.GetSequence(1)) + model = Renumber(aln.GetSequence(1)).CreateFullView() ost.LogInfo(("Computing lDDT between model chain %s and " "reference chain %s") % (ch_mdl, ch_ref)) lddt_scorer = lDDTScorer( references=[reference], model=model, - settings=lddt_settings, - stereochemical_params=stereochemical_parameters) + settings=lddt_settings) try: lddt_results["single_chain_lddt"].append({ "status": "SUCCESS", @@ -529,8 +546,7 @@ def _Main(): qs_scorer.qs_ent_2.ent, qs_scorer.alignments, qs_scorer.calpha_only, - lddt_settings, - stereochemical_parameters) + lddt_settings) lddt_results["oligo_lddt"] = { "status": "SUCCESS", "error": "", diff --git a/modules/mol/alg/pymod/qsscoring.py b/modules/mol/alg/pymod/qsscoring.py index c57ff0820912c139121aedcd2bb8d2f987d65461..094d8b3aa701e6554d9a732afa91f2b7d883c7fc 100644 --- a/modules/mol/alg/pymod/qsscoring.py +++ b/modules/mol/alg/pymod/qsscoring.py @@ -828,7 +828,7 @@ def GetContacts(entity, calpha_only, dist_thr=12.0): class OligoLDDTScorer(object): """A simple class to calculate oligomeric lDDT score.""" - def __init__(self, ref, mdl, alignments, calpha_only, settings, stereochemical_params): + def __init__(self, ref, mdl, alignments, calpha_only, settings): if mdl.chain_count > ref.chain_count: LogWarning('MODEL contains more chains than REFERENCE, ' 'lDDT is not considering them') @@ -838,15 +838,13 @@ class OligoLDDTScorer(object): self.calpha_only = calpha_only self.alignments = alignments self.settings = settings - self.stereochemical_params = stereochemical_params self._lddt = None self._lddt_ref = None self._lddt_mdl = None self.scorer = lDDTScorer( - references=[self.lddt_ref], - model=self.lddt_mdl, - settings=self.settings, - stereochemical_params=self.stereochemical_params) + references=[self.lddt_ref.Select("")], + model=self.lddt_mdl.Select(""), + settings=self.settings) @property def lddt_ref(self): diff --git a/modules/mol/alg/pymod/wrap_mol_alg.cc b/modules/mol/alg/pymod/wrap_mol_alg.cc index 6ec8be1c33a46af1b8e1353476e1e399b9778fc7..efa15a5ea5c7302050d5612c8f7efed024c7c02b 100644 --- a/modules/mol/alg/pymod/wrap_mol_alg.cc +++ b/modules/mol/alg/pymod/wrap_mol_alg.cc @@ -211,21 +211,21 @@ object lDDTScorerInitWrapper(tuple args, dict kwargs){ object self = args[0]; args = tuple(args.slice(1, len(args))); - std::vector<ost::mol::EntityHandle> reference_list_vector; + std::vector<ost::mol::EntityView> reference_list_vector; if(kwargs.contains("references")){ list reference_list = boost::python::extract<list>(kwargs["references"]); int reference_list_length = boost::python::extract<int>(reference_list.attr("__len__")()); for (int i=0; i<reference_list_length; i++) { - reference_list_vector.push_back(boost::python::extract<ost::mol::EntityHandle>(reference_list[i])); + reference_list_vector.push_back(boost::python::extract<ost::mol::EntityView>(reference_list[i])); } kwargs["references"].del(); } else { throw std::invalid_argument("'references' argument is required"); } - ost::mol::EntityHandle model; + ost::mol::EntityView model; if(kwargs.contains("model")){ - model = boost::python::extract<ost::mol::EntityHandle>(kwargs["model"]); + model = boost::python::extract<ost::mol::EntityView>(kwargs["model"]); kwargs["model"].del(); } else { throw std::invalid_argument("'model' argument is required"); @@ -294,7 +294,7 @@ ost::mol::alg::GlobalRDMap prepare_lddt_global_rdmap_wrapper(const list& referen return mol::alg::PreparelDDTGlobalRDMap(reference_list_vector, cutoff_list_vector, sequence_separation, max_dist); } -list get_lddt_per_residue_stats_wrapper(mol::EntityHandle& model, +list get_lddt_per_residue_stats_wrapper(mol::EntityView& model, ost::mol::alg::GlobalRDMap& distance_list, bool structural_checks, String label) { @@ -424,7 +424,7 @@ BOOST_PYTHON_MODULE(_ost_mol_alg) class_<mol::alg::lDDTScorer>("lDDTScorer", no_init) .def("__init__", raw_function(lDDTScorerInitWrapper)) - .def(init<std::vector<mol::EntityHandle>&, mol::EntityHandle&, mol::alg::lDDTSettings&, mol::alg::StereoChemicalProps&>()) + .def(init<std::vector<mol::EntityView>&, mol::EntityView&, mol::alg::lDDTSettings&, mol::alg::StereoChemicalProps&>()) .add_property("global_score", &mol::alg::lDDTScorer::GetGlobalScore) .add_property("conserved_contacts", &mol::alg::lDDTScorer::GetNumConservedContacts) .add_property("total_contacts", &mol::alg::lDDTScorer::GetNumTotalContacts) diff --git a/modules/mol/alg/src/lddt.cc b/modules/mol/alg/src/lddt.cc index e9117757b0696ca982561ca8f82dd5aad328f491..29d1aabfb395304464c8013a2207c590dd7902aa 100644 --- a/modules/mol/alg/src/lddt.cc +++ b/modules/mol/alg/src/lddt.cc @@ -356,7 +356,8 @@ int main (int argc, char **argv) // prints the residue-by-residue statistics std::vector<lDDTLocalScore> local_scores; - local_scores = GetlDDTPerResidueStats(model, glob_dist_list, settings.structural_checks, settings.label); + EntityView outview = model.GetChainList()[0].Select("peptide=true"); + local_scores = GetlDDTPerResidueStats(outview, glob_dist_list, settings.structural_checks, settings.label); PrintlDDTPerResidueStats(local_scores, settings.structural_checks, settings.cutoffs.size()); } return 0; diff --git a/modules/mol/alg/src/local_dist_diff_test.cc b/modules/mol/alg/src/local_dist_diff_test.cc index ff745714d213692be233e011319cc1ebb3908c6c..676c955e10651a7b7bc70165b9da24d1c133e095 100644 --- a/modules/mol/alg/src/local_dist_diff_test.cc +++ b/modules/mol/alg/src/local_dist_diff_test.cc @@ -517,81 +517,56 @@ String lDDTLocalScore::GetHeader(bool structural_checks, int cutoffs_length) { return outstr.str(); } - -lDDTScorer::lDDTScorer(std::vector<EntityHandle>& init_references, - ost::mol::EntityHandle& init_model, +lDDTScorer::lDDTScorer(std::vector<EntityView>& init_references, + ost::mol::EntityView& init_model, lDDTSettings& init_settings): - references(init_references), - model(init_model), settings(init_settings), - _score_calculated(false), - _score_valid(false), - _has_local_scores(false), - _num_cons_con(-1), - _num_tot_con(-1), - _global_score(-1.0) + references_view(init_references), + model_view(init_model) { - model_view = model.GetChainList()[0].Select("peptide=true"); - _PrepareReferences(); - _PrepareGlobalRDMap(); - if (settings.structural_checks) { - if (!stereochemical_params.is_valid) { - throw std::runtime_error( - "Structural checks are ON. Please provide stereochemical_params to enable structural checks."); - } - try { - CheckStructure(model_view, - stereochemical_params.bond_table, - stereochemical_params.angle_table, - stereochemical_params.nonbonded_table, - settings.bond_tolerance, - settings.angle_tolerance); - } catch (std::exception& e) { - std::cout << e.what() << std::endl; - std::stringstream errstr; - errstr << "Structure check failed: " << e.what(); - throw std::runtime_error(errstr.str()); - } - } + _Init(); } -lDDTScorer::lDDTScorer(std::vector<EntityHandle>& init_references, - ost::mol::EntityHandle& init_model, +lDDTScorer::lDDTScorer(std::vector<EntityView>& init_references, + ost::mol::EntityView& init_model, lDDTSettings& init_settings, StereoChemicalProps& init_stereochemical_params): - references(init_references), - model(init_model), settings(init_settings), - stereochemical_params(init_stereochemical_params), - _score_calculated(false), - _score_valid(false), - _has_local_scores(false), - _num_cons_con(-1), - _num_tot_con(-1), - _global_score(-1.0) { - model_view = model.GetChainList()[0].Select("peptide=true"); - _PrepareReferences(); - _PrepareGlobalRDMap(); - if (settings.structural_checks) { - if (!stereochemical_params.is_valid) { - throw std::runtime_error( - "Please provide stereochemical_params to enable structural_checks."); - } - try { - CheckStructure(model_view, - stereochemical_params.bond_table, - stereochemical_params.angle_table, - stereochemical_params.nonbonded_table, - settings.bond_tolerance, - settings.angle_tolerance); - } catch (std::exception& e) { - std::cout << e.what() << std::endl; - std::stringstream errstr; - errstr << "Structure check failed: " << e.what(); - throw std::runtime_error(errstr.str()); - } - } + references_view(init_references), + model_view(init_model), + stereochemical_params(init_stereochemical_params) { + _Init(); + } + +void lDDTScorer::_Init(){ + _score_calculated = false; + _score_valid = false; + _has_local_scores = false; + _num_cons_con = -1; + _num_tot_con = -1; + _global_score = -1.0; + CleanlDDTReferences(references_view); + _PrepareGlobalRDMap(); + if (settings.structural_checks) { + if (!stereochemical_params.is_valid) { + throw std::runtime_error( + "Please provide stereochemical_params to enable structural_checks."); + } + try { + CheckStructure(model_view, + stereochemical_params.bond_table, + stereochemical_params.angle_table, + stereochemical_params.nonbonded_table, + settings.bond_tolerance, + settings.angle_tolerance); + } catch (std::exception& e) { + std::cout << e.what() << std::endl; + std::stringstream errstr; + errstr << "Structure check failed: " << e.what(); + throw std::runtime_error(errstr.str()); + } } +} Real lDDTScorer::GetGlobalScore(){ if (!_score_calculated) { @@ -630,7 +605,7 @@ void lDDTScorer::PrintPerResidueStats(){ settings.cutoffs.size()); } -void lDDTScorer::_PrepareReferences(){ +void lDDTScorer::_PrepareReferences(std::vector<EntityHandle>& references){ for (unsigned int i = 0; i < references.size(); i++) { if (settings.sel != ""){ std::cout << "Performing \"" << settings.sel << "\" selection on reference " << references[i].GetName() << std::endl; @@ -646,8 +621,6 @@ void lDDTScorer::_PrepareReferences(){ references_view.push_back(references[i].CreateFullView()); } } - - CleanlDDTReferences(references_view); } void lDDTScorer::_PrepareGlobalRDMap(){ @@ -708,7 +681,7 @@ void lDDTScorer::_GetLocallDDT(){ if (!_score_calculated){ _ComputelDDT(); } - _local_scores = GetlDDTPerResidueStats(model, + _local_scores = GetlDDTPerResidueStats(model_view, glob_dist_list, settings.structural_checks, settings.label); @@ -1022,12 +995,12 @@ void CheckStructure(EntityView& ent, std::cout << "Distances shorter than tolerance are on average shorter by: " << std::fixed << std::setprecision(5) << clash_info.GetAverageOffset() << std::endl; } -std::vector<lDDTLocalScore> GetlDDTPerResidueStats(EntityHandle& model, +std::vector<lDDTLocalScore> GetlDDTPerResidueStats(EntityView& model, GlobalRDMap& glob_dist_list, bool structural_checks, String label){ std::vector<lDDTLocalScore> scores; - EntityView outv = model.GetChainList()[0].Select("peptide=true"); + EntityView outv = model; for (ChainViewList::const_iterator ci = outv.GetChainList().begin(), ce = outv.GetChainList().end(); ci != ce; ++ci) { for (ResidueViewList::const_iterator rit = ci->GetResidueList().begin(), diff --git a/modules/mol/alg/src/local_dist_diff_test.hh b/modules/mol/alg/src/local_dist_diff_test.hh index cf77b6c495e9afab94dc4c3bcecd044663bedeba..e155f749cacad95353d7362164b71c5b245b3317 100644 --- a/modules/mol/alg/src/local_dist_diff_test.hh +++ b/modules/mol/alg/src/local_dist_diff_test.hh @@ -96,18 +96,16 @@ class lDDTScorer { public: lDDTSettings settings; - ost::mol::EntityHandle model; - std::vector<EntityHandle> references; EntityView model_view; std::vector<EntityView> references_view; GlobalRDMap glob_dist_list; StereoChemicalProps stereochemical_params; - lDDTScorer(std::vector<EntityHandle>& init_references, - ost::mol::EntityHandle& init_model, + lDDTScorer(std::vector<EntityView>& init_references, + ost::mol::EntityView& init_model, lDDTSettings& init_settings); - lDDTScorer(std::vector<EntityHandle>& init_references, - ost::mol::EntityHandle& init_model, + lDDTScorer(std::vector<EntityView>& init_references, + ost::mol::EntityView& init_model, lDDTSettings& init_settings, StereoChemicalProps& init_stereochemical_params); Real GetGlobalScore(); @@ -127,10 +125,11 @@ class lDDTScorer int _num_tot_con; Real _global_score; std::vector<lDDTLocalScore> _local_scores; + void _Init(); void _ComputelDDT(); void _GetLocallDDT(); void _CheckConsistency(); - void _PrepareReferences(); + void _PrepareReferences(std::vector<EntityHandle>& references); void _PrepareGlobalRDMap(); }; @@ -270,7 +269,7 @@ void DLLEXPORT_OST_MOL_ALG CheckStructure(EntityView& ent, Real bond_tolerance, Real angle_tolerance); -std::vector<lDDTLocalScore> DLLEXPORT_OST_MOL_ALG GetlDDTPerResidueStats(EntityHandle& model, +std::vector<lDDTLocalScore> DLLEXPORT_OST_MOL_ALG GetlDDTPerResidueStats(EntityView& model, GlobalRDMap& glob_dist_list, bool structural_checks, String label);