From 92daf6c389b9771a0274f869bacef13d0efdd09e Mon Sep 17 00:00:00 2001
From: Niklaus Johner <nij2003@med.cornell.edu>
Date: Thu, 11 Aug 2011 13:42:08 -0400
Subject: [PATCH] Added a function to evaluate the helicity of a protein
segment on the CoorfFrame and CoordGroup level.
---
.../alg/pymod/export_trajectory_analysis.cc | 1 +
modules/mol/alg/src/trajectory_analysis.cc | 18 +++
modules/mol/alg/src/trajectory_analysis.hh | 1 +
modules/mol/base/pymod/export_coord_frame.cc | 3 +-
modules/mol/base/src/coord_frame.cc | 111 ++++++++++++++++++
modules/mol/base/src/coord_frame.hh | 5 +
6 files changed, 138 insertions(+), 1 deletion(-)
diff --git a/modules/mol/alg/pymod/export_trajectory_analysis.cc b/modules/mol/alg/pymod/export_trajectory_analysis.cc
index d53559215..b2c9dcd56 100644
--- a/modules/mol/alg/pymod/export_trajectory_analysis.cc
+++ b/modules/mol/alg/pymod/export_trajectory_analysis.cc
@@ -48,4 +48,5 @@ void export_TrajectoryAnalysis()
def("AnalyzeAromaticRingInteraction", &AnalyzeAromaticRingInteraction, (arg("traj"), arg("view_ring1"), arg("view_ring2"), arg("stride")=1));
def("AnalyzeAlphaHelixAxis", &AnalyzeAlphaHelixAxis, (arg("traj"), arg("protein_segment"), arg("directions"), arg("centers"), arg("stride")=1));
def("AnalyzeBestFitLine", &AnalyzeBestFitLine, (arg("traj"), arg("protein_segment"), arg("directions"), arg("centers"), arg("stride")=1));
+ def("AnalyzeHelicity", &AnalyzeHelicity, (arg("traj"), arg("protein_segment"), arg("stride")=1));
}
diff --git a/modules/mol/alg/src/trajectory_analysis.cc b/modules/mol/alg/src/trajectory_analysis.cc
index 89b46a525..1688eed33 100644
--- a/modules/mol/alg/src/trajectory_analysis.cc
+++ b/modules/mol/alg/src/trajectory_analysis.cc
@@ -275,6 +275,24 @@ std::vector<Real> AnalyzeAromaticRingInteraction(const CoordGroupHandle& traj, c
return;
}
+ std::vector<Real> AnalyzeHelicity(const CoordGroupHandle& traj, const EntityView& prot_seg,
+ unsigned int stride)
+ {
+ CheckHandleValidity(traj);
+ if (prot_seg.GetAtomCount()==0){
+ throw std::runtime_error("EntityView is empty");
+ }
+ std::vector<unsigned long> indices_c,indices_o, indices_n, indices_ca;
+ std::vector<Real> helicity;
+ helicity.reserve(ceil(traj.GetFrameCount()/float(stride)));
+ GetCaCONIndices(prot_seg, indices_ca, indices_c, indices_o, indices_n);
+ for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
+ CoordFramePtr frame=traj.GetFrame(i);
+ helicity.push_back(frame->GetAlphaHelixContent(indices_ca,indices_c,indices_o,indices_n));
+ }
+ return helicity;
+ }
+
}}} //ns
diff --git a/modules/mol/alg/src/trajectory_analysis.hh b/modules/mol/alg/src/trajectory_analysis.hh
index 6837e9fb3..5c2ec6de2 100644
--- a/modules/mol/alg/src/trajectory_analysis.hh
+++ b/modules/mol/alg/src/trajectory_analysis.hh
@@ -45,5 +45,6 @@ namespace ost { namespace mol { namespace alg {
std::vector<Real> DLLEXPORT_OST_MOL_ALG AnalyzeAromaticRingInteraction(const CoordGroupHandle& traj, const EntityView& view_ring1, const EntityView& view_ring2,unsigned int stride=1);
void DLLEXPORT_OST_MOL_ALG AnalyzeAlphaHelixAxis(const CoordGroupHandle& traj, const EntityView& prot_seg, geom::Vec3List& directions, geom::Vec3List& centers, unsigned int stride=1);
void DLLEXPORT_OST_MOL_ALG AnalyzeBestFitLine(const CoordGroupHandle& traj, const EntityView& prot_seg, geom::Vec3List& directions, geom::Vec3List& centers, unsigned int stride=1);
+ std::vector<Real> DLLEXPORT_OST_MOL_ALG AnalyzeHelicity(const CoordGroupHandle& traj, const EntityView& prot_seg, unsigned int stride=1);
}}}//ns
#endif
diff --git a/modules/mol/base/pymod/export_coord_frame.cc b/modules/mol/base/pymod/export_coord_frame.cc
index 5ec6f157b..6636421ad 100644
--- a/modules/mol/base/pymod/export_coord_frame.cc
+++ b/modules/mol/base/pymod/export_coord_frame.cc
@@ -34,7 +34,7 @@ geom::Vec3 (CoordFrame::*get_cm)(const mol::EntityView& sele) = &CoordFrame::Get
Real (CoordFrame::*get_dist_cm)(const mol::EntityView& sele1, const mol::EntityView& sele2) = &CoordFrame::GetDistanceBetwCenterOfMass;
Real (CoordFrame::*get_rmsd)(const mol::EntityView& Reference_View, const mol::EntityView& sele_View) = &CoordFrame::GetRMSD;
Real (CoordFrame::*get_min_dist)(const mol::EntityView& view1, const mol::EntityView& view2) = &CoordFrame::GetMinDistance;
-
+Real (CoordFrame::*get_alpha)(const mol::EntityView& segment) = &CoordFrame::GetAlphaHelixContent;
void export_CoordFrame()
{
class_<CoordFrame>("CoordFrame",no_init)
@@ -47,6 +47,7 @@ void export_CoordFrame()
.def("GetRMSD",get_rmsd)
.def("GetMinDistance",get_min_dist)
.def("GetODRLine",&geom::Vec3List::GetODRLine)
+ .def("GetAlphaHelixContent",get_alpha)
;
def("CreateCoordFrame",CreateCoordFrame);
}
diff --git a/modules/mol/base/src/coord_frame.cc b/modules/mol/base/src/coord_frame.cc
index 01db3f0a2..efeb37e20 100644
--- a/modules/mol/base/src/coord_frame.cc
+++ b/modules/mol/base/src/coord_frame.cc
@@ -28,6 +28,7 @@
#include <ost/mol/view_op.hh>
#include <ost/mol/mol.hh>
#include "coord_frame.hh"
+#include <ost/base.hh>
namespace ost { namespace mol {
@@ -101,6 +102,35 @@ namespace ost { namespace mol {
}
}
+ void GetCaCONIndices(const EntityView& segment, std::vector<unsigned long>& indices_ca,
+ std::vector<unsigned long>& indices_c, std::vector<unsigned long>& indices_o, std::vector<unsigned long>& indices_n){
+ mol::AtomView a,n,c,o;
+ ResidueViewList residues=segment.GetResidueList();
+ indices_ca.reserve(residues.size());
+ indices_n.reserve(residues.size());
+ indices_c.reserve(residues.size());
+ indices_o.reserve(residues.size());
+ for (ResidueViewList::const_iterator res=residues.begin()+1,
+ e=residues.end(); res!=e-1; ++res){
+ if (!InSequence((*res).GetHandle(),(*(res+1)).GetHandle())) throw std::runtime_error("Residues are not in a continuous sequence");
+ }
+ for (ResidueViewList::const_iterator res=residues.begin(),
+ e=residues.end(); res!=e; ++res) {
+ a=res->FindAtom("CA");
+ CheckHandleValidity(a);
+ indices_ca.push_back(a.GetIndex());
+ c=res->FindAtom("C");
+ CheckHandleValidity(c);
+ indices_c.push_back(c.GetIndex());
+ o=res->FindAtom("O");
+ CheckHandleValidity(o);
+ indices_o.push_back(o.GetIndex());
+ n=res->FindAtom("N");
+ CheckHandleValidity(n);
+ indices_n.push_back(n.GetIndex());
+ }
+ }
+
void GetMasses(const EntityView& sele, std::vector<Real>& masses){
// This function returns a vector containing the atom masses of the atoms in an EntityView
// It is used together with GetIndices to accelerate the extraction of RMSD from a trajectory
@@ -255,6 +285,87 @@ namespace ost { namespace mol {
return atoms_pos_list.FitCylinder();
}
+ Real CoordFrame::GetAlphaHelixContent(const mol::EntityView& segment){
+ //Returns the percentage of residues in the EntityView (segment) that are in an alpha-helix
+ //Each residue is assigned as being in an alpha-helix or not, based on and backbone torsion angles
+ //and the hydrogen bonding (O-N) (this term is used only if the segment contains at least 8 residues
+ //The entity view should be a single segment with no gaps and no missing N,CA,C,O backbone atoms
+ //the first and last residues will not be asessed for helicity
+ CheckHandleValidity(segment);
+ std::vector<unsigned long> indices_c, indices_n, indices_ca, indices_o;
+ GetCaCONIndices(segment, indices_ca, indices_c, indices_o, indices_n);
+ return CoordFrame::GetAlphaHelixContent(indices_ca,indices_c,indices_o,indices_n);
+ }
+
+ Real Eval2AngleDist(Real phi, Real phi0, Real psi, Real psi0, Real delta){
+ //Returns a score between 0 and 1 measuring the distance between
+ //a pair of angles and a pair of reference angles
+ Real d1,d2,d;
+ d1=abs(phi-phi0);
+ d2=abs(psi-psi0);
+ if (d1>M_PI) d1=abs(d1-2.*M_PI);
+ if (d2>M_PI) d1=abs(d2-2.*M_PI);
+ d=(pow(d1,2.)+pow(d2,2.))/pow(delta,2.);
+ return 1./(1+d);
+ }
+
+ Real CoordFrame::GetAlphaHelixContent(std::vector<unsigned long>& indices_ca, std::vector<unsigned long>& indices_c,
+ std::vector<unsigned long>& indices_o, std::vector<unsigned long>& indices_n){
+ //See CoordFrame::GetAlphaHelixContent(const mol::EntityView) above.
+ unsigned long n_atoms=indices_ca.size();
+ geom::Vec3 c_previous,n,ca,c,n_next;
+ Real phi,psi,phi_0=-1.2,psi_0=-0.785,delta=0.8,d_0=3.3;
+ unsigned long n_helical_res=0;
+ std::vector<Real> score,dist,score2,dist2;
+ score.reserve(n_atoms-2);
+ score2.reserve(n_atoms-2);
+ dist.reserve(n_atoms-2);
+ dist2.reserve(n_atoms-2);
+ if (n_atoms!=indices_n.size()||n_atoms!=indices_c.size()||n_atoms!=indices_o.size()){
+ throw std::runtime_error("not same numbers of CA, C, O and N atoms in the selection");
+ }
+ if (n_atoms<=5){
+ throw std::runtime_error("At least five residues are needed to calulate an alpha helix similarity");
+ }
+ c=(*this)[indices_c[0]];
+ n_next=(*this)[indices_n[1]];
+ for (unsigned long i=1; i!=n_atoms-1; ++i){
+ c_previous=c;
+ n=n_next;
+ ca=(*this)[indices_ca[i]];
+ c=(*this)[indices_c[i]];
+ n_next=(*this)[indices_n[i+1]];
+ phi=geom::DihedralAngle(c_previous,n,ca,c);
+ psi=geom::DihedralAngle(n,ca,c,n_next);
+ score.push_back(Eval2AngleDist(phi,phi_0,psi,psi_0,delta));
+ }
+ score2[0]=pow(score[0]*score[1],3./2.);
+ score2[n_atoms-3]=pow(score[n_atoms-3]*score[n_atoms-4],3./2.);
+ for (unsigned long i=1; i!=n_atoms-3; ++i){
+ score2[i]=score[i-1]*score[i]*score[i+1];
+ }
+ if (n_atoms>=8){
+ for (unsigned long i=1; i!=n_atoms-1; ++i){
+ if (i<n_atoms-4){
+ dist.push_back(geom::Distance((*this)[indices_o[i]],(*this)[indices_n[i+4]]));
+ if (i>=5) {
+ dist2.push_back(std::min(dist[i-1],dist[i-5]));
+ }
+ else dist2.push_back(dist[i-1]);
+ }
+ else dist2.push_back(geom::Distance((*this)[indices_n[i]],(*this)[indices_o[i-4]]));
+ }
+ for (unsigned long i=0; i!=n_atoms-2; ++i){
+ if (dist2[i]<=d_0 || score2[i]>=0.3)n_helical_res+=1;
+ }
+ }
+ else {
+ for (unsigned long i=0; i!=n_atoms-2; ++i){
+ if (score2[i]>=0.3)n_helical_res+=1;
+ }
+ }
+ return Real(n_helical_res)/Real(n_atoms-2);
+ }
diff --git a/modules/mol/base/src/coord_frame.hh b/modules/mol/base/src/coord_frame.hh
index fe94ecc9f..6623756f5 100644
--- a/modules/mol/base/src/coord_frame.hh
+++ b/modules/mol/base/src/coord_frame.hh
@@ -62,6 +62,9 @@ public:
Real GetMinDistBetwCenterOfMassAndView(const mol::EntityView& view_cm, const mol::EntityView& view_atoms);
geom::Line3 GetODRLine(std::vector<unsigned long>& indices_ca);
geom::Line3 FitCylinder(std::vector<unsigned long>& indices_ca);
+ Real GetAlphaHelixContent(std::vector<unsigned long>& indices_ca, std::vector<unsigned long>& indices_ca,
+ std::vector<unsigned long>& indices_o, std::vector<unsigned long>& indices_n);
+ Real GetAlphaHelixContent(const mol::EntityView& segment);
};
void GetIndices(const EntityView& sele, std::vector<unsigned long>& indices);
@@ -69,6 +72,8 @@ public:
void GetIndicesAndMasses(const EntityView& sele, std::vector<unsigned long>& indices,std::vector<Real>& masses);
void GetPositions(const EntityView& sele, std::vector<geom::Vec3>& ref_pos);
void GetCaIndices(const EntityView& segment, std::vector<unsigned long>& indices_ca);
+ void GetCaCONIndices(const EntityView& segment, std::vector<unsigned long>& indices_ca, std::vector<unsigned long>& indices_c,
+ std::vector<unsigned long>& indices_o, std::vector<unsigned long>& indices_n);
typedef boost::shared_ptr<CoordFrame> CoordFramePtr;
typedef std::vector<CoordFramePtr> CoordFrameList;
--
GitLab