diff --git a/modules/seq/base/pymod/export_sequence.cc b/modules/seq/base/pymod/export_sequence.cc index ab36fdaed36c6fd40b3d8669b0a135b0a0cef3bd..b3632fc507db5bdf281a80c7384173a70fe7e0d0 100644 --- a/modules/seq/base/pymod/export_sequence.cc +++ b/modules/seq/base/pymod/export_sequence.cc @@ -301,6 +301,7 @@ void export_sequence() .def("GetLength", &AlignmentHandle::GetLength) .def("__len__", &AlignmentHandle::GetLength) .def("GetSequences", &AlignmentHandle::GetSequences) + .def("GetCoverage", &AlignmentHandle::GetSequences) .def("AttachView", attach_view_a) .def("AttachView", attach_view_b) .def("Cut", &AlignmentHandle::Cut) diff --git a/modules/seq/base/src/alignment_handle.cc b/modules/seq/base/src/alignment_handle.cc index d2a7ab85d387051d43edbadac9419925b5ac7fee..bf2c12fdf61f3a16e340891915011011f8c25fa0 100644 --- a/modules/seq/base/src/alignment_handle.cc +++ b/modules/seq/base/src/alignment_handle.cc @@ -263,4 +263,12 @@ int AlignmentHandle::GetSequenceOffset(int seq_index) this->CheckValidity(); return impl_->GetSequence(seq_index)->GetOffset(); } + +Real AlignmentHandle::GetCoverage(int seq_index) const +{ + this->CheckValidity(); + return impl_->GetCoverage(seq_index); +} + + }} diff --git a/modules/seq/base/src/alignment_handle.hh b/modules/seq/base/src/alignment_handle.hh index bcba027d92f718fb0d874e96361dfc4f0387ffbd..823e78be1ef5e2e7ade4e5c0f690df9d1190f261 100644 --- a/modules/seq/base/src/alignment_handle.hh +++ b/modules/seq/base/src/alignment_handle.hh @@ -158,6 +158,13 @@ public: iterator end() const; bool IsValid() const { return impl_.get()!=0; } + + ///\brief get coverage of a specifi sequence + /// + /// returns a value representing how extensively the specified sequence + /// covers the first sequence (sequence 0). The function return a value + /// between 0 (no coverage) and 1 (full coverage) + Real GetCoverage(int seq_index) const; private: void CheckValidity() const; diff --git a/modules/seq/base/src/impl/sequence_list_impl.cc b/modules/seq/base/src/impl/sequence_list_impl.cc index b955d8b4cf1244e144d9c3bed82b4688b47afbed..16759dbcf67dbc1d71d356dfca886b2922dd5fd3 100644 --- a/modules/seq/base/src/impl/sequence_list_impl.cc +++ b/modules/seq/base/src/impl/sequence_list_impl.cc @@ -178,4 +178,22 @@ SequenceListImplPtr SequenceListImpl::Slice(int first, int n) const } } +Real SequenceListImpl::GetCoverage(int seq_index) const +{ + int a=0, b=0; + String seq_string_0 = GetSequence(0)->GetString(); + String seq_string_index = GetSequence(seq_index)->GetString(); + for (int i=0; i<GetSequence(0)->GetLength(); ++i) { + if (seq_string_0[i]!='-') { + a+=1; + if (seq_string_index[i]!='-') { + b+=1; + } + } + } + Real coverage=Real(b)/Real(a); + return coverage; +} + + }}} diff --git a/modules/seq/base/src/impl/sequence_list_impl.hh b/modules/seq/base/src/impl/sequence_list_impl.hh index 66bedd61f47c4cb5d40374b92dfc94c2213ecee9..ab382b8c8faf98687757f7c916c2c92ccc36207a 100644 --- a/modules/seq/base/src/impl/sequence_list_impl.hh +++ b/modules/seq/base/src/impl/sequence_list_impl.hh @@ -78,7 +78,7 @@ public: return this->GetMaxLength()==this->GetMinLength(); } - + Real GetCoverage(int seq_index) const; SequenceListImplPtr Copy() const;