From 67122394ca894dc1dbbdfa731d5c5ceaf32fa68f Mon Sep 17 00:00:00 2001
From: Marco Biasini <marco.biasini@unibas.ch>
Date: Tue, 29 Mar 2011 11:21:42 +0200
Subject: [PATCH] add GetMatchingBackboneViews to AlignmentHandle
---
modules/mol/base/pymod/export_entity_view.cc | 4 ++
modules/mol/base/src/entity_view.hh | 2 +
modules/seq/base/doc/seq.rst | 19 ++++++++++
.../base/pymod/const_seq_list_export_def.hh | 1 +
modules/seq/base/pymod/export_sequence.cc | 7 ++--
modules/seq/base/src/alignment_handle.cc | 37 +++++++++++++++++++
modules/seq/base/src/alignment_handle.hh | 2 +
modules/seq/base/src/impl/sequence_impl.hh | 11 +++++-
modules/seq/base/src/sequence_list.cc | 12 ++++++
modules/seq/base/src/sequence_list.hh | 3 ++
10 files changed, 93 insertions(+), 5 deletions(-)
diff --git a/modules/mol/base/pymod/export_entity_view.cc b/modules/mol/base/pymod/export_entity_view.cc
index 5e3dc12d2..55d408365 100644
--- a/modules/mol/base/pymod/export_entity_view.cc
+++ b/modules/mol/base/pymod/export_entity_view.cc
@@ -20,6 +20,8 @@
#include <boost/python/suite/indexing/vector_indexing_suite.hpp>
using namespace boost::python;
+#include <ost/export_helper/pair_to_tuple_conv.hh>
+
#include <ost/mol/entity_view.hh>
#include <ost/mol/query.hh>
#include <ost/mol/mol.hh>
@@ -184,6 +186,8 @@ void export_EntityView()
.def("GetBounds", &EntityView::GetBounds)
.add_property("bounds", &EntityView::GetBounds)
;
+ to_python_converter<std::pair<EntityView, EntityView>,
+ PairToTupleConverter<EntityView, EntityView> >();
def("Union", &Union);
def("Difference", &Difference);
def("Intersection", &Intersection);
diff --git a/modules/mol/base/src/entity_view.hh b/modules/mol/base/src/entity_view.hh
index 5d2a0848c..6e5b23716 100644
--- a/modules/mol/base/src/entity_view.hh
+++ b/modules/mol/base/src/entity_view.hh
@@ -335,6 +335,8 @@ private:
EntityViewDataPtr data_;
};
+typedef std::pair<EntityView, EntityView> EntityViewPair;
+
}} // ns
#endif // OST_ENTITY_VIEW_HH
diff --git a/modules/seq/base/doc/seq.rst b/modules/seq/base/doc/seq.rst
index ade28f583..f42a2f574 100644
--- a/modules/seq/base/doc/seq.rst
+++ b/modules/seq/base/doc/seq.rst
@@ -357,3 +357,22 @@ an alignment:
If master is set to -1, all sequences in the region are affected, otherwise
only the sequence at index equal to master is shifted.
+
+ .. method:: GetMatchingBackboneViews(index1=0, index2=1)
+
+ Returns a tuple of entity views containing matching backbone atoms for the
+ two sequences at index1 and index2, respectively. For each aligned column in
+ the alignment, backbone atoms are added to the view if both aligned residues
+ have them. It is guaranteed that the two views contain the same number of
+ atoms and that the order of the atoms in the two views is the same.
+
+ The output of this function can be used to superpose two structures with
+ :func:`~ost.mol.alg.SuperposeSVD`.
+
+
+ :param index1: The index of the first sequence
+
+ :param index2: The index of the second sequence.
+
+ :raises: In case one of the two sequences doesn't have an attached view, a
+ :exc:`RuntimeError` is raised.
diff --git a/modules/seq/base/pymod/const_seq_list_export_def.hh b/modules/seq/base/pymod/const_seq_list_export_def.hh
index 159857c3b..cd1f52853 100644
--- a/modules/seq/base/pymod/const_seq_list_export_def.hh
+++ b/modules/seq/base/pymod/const_seq_list_export_def.hh
@@ -26,6 +26,7 @@
.def("IsValid", &C::IsValid) \
.def("Take", &C::Take) \
.def("Slice", &C::Slice) \
+ .def("FindSequence", &C::FindSequence) \
.def("SequencesHaveEqualLength", \
&C::SequencesHaveEqualLength) \
.def("__getitem__", &C::operator[]) \
diff --git a/modules/seq/base/pymod/export_sequence.cc b/modules/seq/base/pymod/export_sequence.cc
index 454245778..c431dfa9b 100644
--- a/modules/seq/base/pymod/export_sequence.cc
+++ b/modules/seq/base/pymod/export_sequence.cc
@@ -21,8 +21,6 @@
#include <boost/python/register_ptr_to_python.hpp>
#include <boost/python/suite/indexing/vector_indexing_suite.hpp>
-
-#include <ost/export_helper/pair_to_tuple_conv.hh>
#include <ost/generic_property.hh>
#include <ost/export_helper/generic_property_def.hh>
#include <ost/info/info.hh>
@@ -299,6 +297,8 @@ void export_sequence()
.def("Copy", &AlignmentHandle::Copy)
.def("ToString", &AlignmentHandle::ToString)
.def("GetLength", &AlignmentHandle::GetLength)
+ .def("GetMatchingBackboneViews", &AlignmentHandle::GetMatchingBackboneViews,
+ (arg("index1")=0, arg("index2")=1))
.def("__len__", &AlignmentHandle::GetLength)
.def("GetSequences", &AlignmentHandle::GetSequences)
.def("GetCoverage", &AlignmentHandle::GetSequences)
@@ -357,8 +357,7 @@ void export_sequence()
.def("__getitem__", &do_slice_b)
;
implicitly_convertible<SequenceList, ConstSequenceList>();
- to_python_converter<std::pair<mol::EntityView, mol::EntityView>,
- PairToTupleConverter<mol::EntityView, mol::EntityView> >();
+
def("CreateSequenceList", &CreateSequenceList);
def("SequenceFromChain", seq_from_chain_a);
def("SequenceFromChain", seq_from_chain_b);
diff --git a/modules/seq/base/src/alignment_handle.cc b/modules/seq/base/src/alignment_handle.cc
index c39856e93..3e4ed4256 100644
--- a/modules/seq/base/src/alignment_handle.cc
+++ b/modules/seq/base/src/alignment_handle.cc
@@ -20,8 +20,12 @@
/*
Author: Marco Biasini
*/
+
#include <ost/invalid_handle.hh>
+#include <ost/mol/residue_view.hh>
+#include <ost/mol/atom_view.hh>
#include <ost/seq/alignment_handle.hh>
+#include <ost/mol/residue_view.hh>
#include <ost/seq/impl/sequence_list_impl.hh>
#include <ost/seq/impl/sequence_impl.hh>
#include <ost/seq/sequence_list.hh>
@@ -270,5 +274,38 @@ Real AlignmentHandle::GetCoverage(int seq_index) const
return impl_->GetCoverage(seq_index);
}
+mol::EntityViewPair AlignmentHandle::GetMatchingBackboneViews(int idx0, int idx1) const
+{
+ this->CheckValidity();
+ const impl::SequenceImpl& s1=*impl_->GetSequence(idx0).get();
+ const impl::SequenceImpl& s2=*impl_->GetSequence(idx1).get();
+ if (!s1.HasAttachedView() || !s2.HasAttachedView()) {
+ throw std::runtime_error("both sequences must have a view attached");
+ }
+ mol::EntityView v1=s1.GetAttachedView().CreateEmptyView();
+ mol::EntityView v2=s2.GetAttachedView().CreateEmptyView();
+ for (int i=0; i<s1.GetLength(); ++i) {
+ if (s1[i]=='-' && s2[i]=='-') {
+ continue;
+ }
+ mol::ResidueView r1=s1.GetResidue(i);
+ mol::ResidueView r2=s2.GetResidue(i);
+ if (!r1.IsValid() || !r2.IsValid()) {
+ continue;
+ }
+ const char* bb_anames[]={"N", "CA", "C", "O"};
+ //for (size_t j=0; )
+ for (size_t j=0; j<4; ++j) {
+ mol::AtomView a1=r1.FindAtom(bb_anames[j]);
+ mol::AtomView a2=r2.FindAtom(bb_anames[j]);
+ if (!a1.IsValid() || !a2.IsValid()) {
+ continue;
+ }
+ v1.AddAtom(a1);
+ v2.AddAtom(a2);
+ }
+ }
+ return mol::EntityViewPair(v1, v2);
+}
}}
diff --git a/modules/seq/base/src/alignment_handle.hh b/modules/seq/base/src/alignment_handle.hh
index 151506d25..9a109215b 100644
--- a/modules/seq/base/src/alignment_handle.hh
+++ b/modules/seq/base/src/alignment_handle.hh
@@ -107,6 +107,8 @@ public:
void AttachView(int seq_index, const mol::EntityView& view,
const String& chain_name);
+
+ mol::EntityViewPair GetMatchingBackboneViews(int idx0=0, int idx1=1) const;
/// \brief set name of sequence
void SetSequenceName(int seq_index, const String& name);
diff --git a/modules/seq/base/src/impl/sequence_impl.hh b/modules/seq/base/src/impl/sequence_impl.hh
index fd7fd0dfe..beae61e03 100644
--- a/modules/seq/base/src/impl/sequence_impl.hh
+++ b/modules/seq/base/src/impl/sequence_impl.hh
@@ -137,6 +137,15 @@ public:
void Append(char olc);
+ char& operator[](size_t index)
+ {
+ return seq_string_[index];
+ }
+ char operator[](size_t index) const
+ {
+ return seq_string_[index];
+ }
+
private:
/// \brief Recalculates gap shifts from sequence.
@@ -160,7 +169,7 @@ private:
std::list<Shift> shifts_;
bool editing_;
int offset_;
- mol::EntityView attached_view_;
+ mol::EntityView attached_view_;
};
/// \internal
diff --git a/modules/seq/base/src/sequence_list.cc b/modules/seq/base/src/sequence_list.cc
index 699a7a0eb..eb45568fd 100644
--- a/modules/seq/base/src/sequence_list.cc
+++ b/modules/seq/base/src/sequence_list.cc
@@ -72,6 +72,12 @@ SequenceList::Iterator SequenceList::End() const
return SequenceList::Iterator(impl_->End(), impl_->End());
}
+SequenceHandle SequenceList::FindSequence(const String& name) const
+{
+ this->CheckValidity();
+ return SequenceHandle(impl_->FindSequence(name));
+}
+
bool SequenceList::IsValid() const
{
return impl_.get()!=NULL;
@@ -226,6 +232,12 @@ ConstSequenceList ConstSequenceList::Slice(int first, int n) const
return ConstSequenceList(impl_->Slice(first, n));
}
+ConstSequenceHandle ConstSequenceList::FindSequence(const String& name) const
+{
+ this->CheckValidity();
+ return ConstSequenceHandle(impl_->FindSequence(name));
+}
+
ConstSequenceList DLLEXPORT_OST_SEQ CreateConstSequenceList()
{
return ConstSequenceList(impl::SequenceListImplPtr(new impl::SequenceListImpl));
diff --git a/modules/seq/base/src/sequence_list.hh b/modules/seq/base/src/sequence_list.hh
index 62eea1041..71067438d 100644
--- a/modules/seq/base/src/sequence_list.hh
+++ b/modules/seq/base/src/sequence_list.hh
@@ -61,6 +61,8 @@ public:
/// \brief create a sequence list from the given splice interval
ConstSequenceList Slice(int first, int n) const;
+
+ ConstSequenceHandle FindSequence(const String& name) const;
int GetMinLength() const;
int GetMaxLength() const;
/// \internal
@@ -107,6 +109,7 @@ public:
int GetMinLength() const;
int GetMaxLength() const;
+ SequenceHandle FindSequence(const String& name) const;
// \internal
impl::SequenceListImplPtr& Impl() const;
SequenceList(const impl::SequenceListImplPtr& impl);
--
GitLab