Something went wrong on our end
-
Niklaus Johner authored
of the fitted cylinder. The input was also changed from taking both a direction and center as initial guess to just the direction. The center is guessed from the geometrical center of the atoms and fitted during the optimization procedure.
Niklaus Johner authoredof the fitted cylinder. The input was also changed from taking both a direction and center as initial guess to just the direction. The center is guessed from the geometrical center of the atoms and fitted during the optimization procedure.
trajectory_analysis.cc 16.29 KiB
//------------------------------------------------------------------------------
// This file is part of the OpenStructure project <www.openstructure.org>
//
// Copyright (C) 2008-2011 by the OpenStructure authors
//
// This library is free software; you can redistribute it and/or modify it under
// the terms of the GNU Lesser General Public License as published by the Free
// Software Foundation; either version 3.0 of the License, or (at your option)
// any later version.
// This library is distributed in the hope that it will be useful, but WITHOUT
// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
// details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with this library; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
//------------------------------------------------------------------------------
/*
* Author Niklaus Johner
*/
#include <stdexcept>
#include <ost/base.hh>
#include <ost/geom/vec3.hh>
#include <ost/base.hh>
#include <ost/geom/geom.hh>
#include <ost/mol/mol.hh>
#include <ost/mol/view_op.hh>
#include "trajectory_analysis.hh"
namespace ost { namespace mol { namespace alg {
geom::Vec3List AnalyzeAtomPos(const CoordGroupHandle& traj, const AtomHandle& a1, unsigned int stride)
//This function extracts the position of an atom from a trajectory and returns it as a vector of geom::Vec3
{
CheckHandleValidity(traj);
geom::Vec3List pos;
pos.reserve(ceil(traj.GetFrameCount()/float(stride)));
int i1=a1.GetIndex();
for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
CoordFramePtr frame=traj.GetFrame(i);
pos.push_back(frame->GetAtomPos(i1));
}
return pos;
}
std::vector<Real> AnalyzeDistanceBetwAtoms(const CoordGroupHandle& traj, const AtomHandle& a1, const AtomHandle& a2,
unsigned int stride)
//This function extracts the distance between two atoms from a trajectory and returns it as a vector
{
CheckHandleValidity(traj);
std::vector<Real> dist;
dist.reserve(ceil(traj.GetFrameCount()/float(stride)));
int i1=a1.GetIndex();
int i2=a2.GetIndex();
for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
CoordFramePtr frame=traj.GetFrame(i);
dist.push_back(frame->GetDistanceBetwAtoms(i1,i2));
}
return dist;
}
std::vector<Real> AnalyzeAngle(const CoordGroupHandle& traj, const AtomHandle& a1, const AtomHandle& a2,
const AtomHandle& a3, unsigned int stride)
//This function extracts the angle between three atoms from a trajectory and returns it as a vector
{
CheckHandleValidity(traj);
std::vector<Real> ang;
ang.reserve(ceil(traj.GetFrameCount()/float(stride)));
int i1=a1.GetIndex(),i2=a2.GetIndex(),i3=a3.GetIndex();
for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
CoordFramePtr frame=traj.GetFrame(i);
ang.push_back(frame->GetAngle(i1,i2,i3));
}
return ang;
}
std::vector<Real> AnalyzeDihedralAngle(const CoordGroupHandle& traj, const AtomHandle& a1, const AtomHandle& a2,
const AtomHandle& a3, const AtomHandle& a4, unsigned int stride)
//This function extracts the dihedral angle between four atoms from a trajectory and returns it as a vector
{
CheckHandleValidity(traj);
std::vector<Real> ang;
ang.reserve(ceil(traj.GetFrameCount()/float(stride)));
int i1=a1.GetIndex(),i2=a2.GetIndex(),i3=a3.GetIndex(),i4=a4.GetIndex();
for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
CoordFramePtr frame=traj.GetFrame(i);
ang.push_back(frame->GetDihedralAngle(i1,i2,i3,i4));
}
return ang;
}
geom::Vec3List AnalyzeCenterOfMassPos(const CoordGroupHandle& traj, const EntityView& sele,unsigned int stride)
//This function extracts the position of the CenterOfMass of a selection (entity view) from a trajectory
//and returns it as a vector.
{
CheckHandleValidity(traj);
geom::Vec3List pos;
pos.reserve(ceil(traj.GetFrameCount()/float(stride)));
std::vector<unsigned long> indices;
std::vector<Real> masses;
GetIndicesAndMasses(sele, indices, masses);
for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
CoordFramePtr frame=traj.GetFrame(i);
pos.push_back(frame->GetCenterOfMassPos(indices,masses));
}
return pos;
}
std::vector<Real> AnalyzeDistanceBetwCenterOfMass(const CoordGroupHandle& traj, const EntityView& sele1,
const EntityView& sele2, unsigned int stride)
//This function extracts the distance between the CenterOfMass of two selection (entity views) from a trajectory
//and returns it as a vector.
{
CheckHandleValidity(traj);
std::vector<Real> dist;
dist.reserve(ceil(traj.GetFrameCount()/float(stride)));
std::vector<unsigned long> indices1,indices2;
std::vector<Real> masses1,masses2;
GetIndicesAndMasses(sele1, indices1, masses1);
GetIndicesAndMasses(sele2, indices2, masses2);
for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
CoordFramePtr frame=traj.GetFrame(i);
dist.push_back(frame->GetDistanceBetwCenterOfMass(indices1,masses1,indices2,masses2));
}
return dist;
}
std::vector<Real> AnalyzeRMSD(const CoordGroupHandle& traj, const EntityView& reference_view,
const EntityView& sele_view, unsigned int stride)
// This function extracts the rmsd between two entity views and returns it as a vector
// The views don't have to be from the same entity
// If you want to compare to frame i of the trajectory t, first use t.CopyFrame(i) for example:
// eh=io.LoadPDB(...),t=io.LoadCHARMMTraj(eh,...);Sele=eh.Select(...);t.CopyFrame(0);mol.alg.AnalyzeRMSD(t,Sele,Sele)
{
CheckHandleValidity(traj);
int count_ref=reference_view.GetAtomCount();
int count_sele=sele_view.GetAtomCount();
if (count_ref!=count_sele){
throw Error("atom counts of the two views are not equal");
}
std::vector<Real> rmsd;
rmsd.reserve(ceil(traj.GetFrameCount()/float(stride)));
std::vector<unsigned long> sele_indices;
std::vector<geom::Vec3> ref_pos;
GetIndices(sele_view, sele_indices);
GetPositions(reference_view, ref_pos);
for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
CoordFramePtr frame=traj.GetFrame(i);
rmsd.push_back(frame->GetRMSD(ref_pos,sele_indices));
}
return rmsd;
}
std::vector<Real> AnalyzeMinDistance(const CoordGroupHandle& traj, const EntityView& view1,
const EntityView& view2,unsigned int stride)
// This function extracts the minimal distance between two sets of atoms (view1 and view2) for
// each frame in a trajectory and returns it as a vector.
{
CheckHandleValidity(traj);
if (!view1.HasAtoms()){
throw Error("first EntityView is empty");
}
if (!view2.HasAtoms()){
throw Error("second EntityView is empty");
}
std::vector<Real> dist;
dist.reserve(ceil(traj.GetFrameCount()/float(stride)));
std::vector<unsigned long> indices1,indices2;
GetIndices(view1, indices1);
GetIndices(view2, indices2);
for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
CoordFramePtr frame=traj.GetFrame(i);
dist.push_back(frame->GetMinDistance(indices1,indices2));
}
return dist;
}
std::vector<Real> AnalyzeMinDistanceBetwCenterOfMassAndView(const CoordGroupHandle& traj, const EntityView& view_cm,
const EntityView& view_atoms,unsigned int stride)
// This function extracts the minimal distance between a set of atoms (view_atoms) and the center of mass
// of a second set of atoms (view_cm) for each frame in a trajectory and returns it as a vector.
{
CheckHandleValidity(traj);
if (!view_cm.HasAtoms()){
throw Error("first EntityView is empty");
}
if (!view_atoms.HasAtoms()){
throw Error("second EntityView is empty");
}
std::vector<Real> dist, masses_cm;
dist.reserve(ceil(traj.GetFrameCount()/float(stride)));
std::vector<unsigned long> indices_cm,indices_atoms;
GetIndicesAndMasses(view_cm, indices_cm,masses_cm);
GetIndices(view_atoms, indices_atoms);
for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
CoordFramePtr frame=traj.GetFrame(i);
dist.push_back(frame->GetMinDistBetwCenterOfMassAndView(indices_cm, masses_cm, indices_atoms));
}
return dist;
}
std::vector<Real> AnalyzeAromaticRingInteraction(const CoordGroupHandle& traj, const EntityView& view_ring1,
const EntityView& view_ring2,unsigned int stride)
// This function is a crude analysis of aromatic ring interactions. For each frame in a trajectory, it calculates
// the minimal distance between the atoms in one view and the center of mass of the other
// and vice versa, and returns the minimum between these two minimal distances.
// Concretely, if the two views are the heavy atoms of two rings, then it returns the minimal
// center of mass - heavy atom distance betweent he two rings
{
CheckHandleValidity(traj);
if (!view_ring1.HasAtoms()){
throw Error("first EntityView is empty");
}
if (!view_ring2.HasAtoms()){
throw Error("second EntityView is empty");
}
std::vector<Real> dist, masses_ring1,masses_ring2;
dist.reserve(ceil(traj.GetFrameCount()/float(stride)));
std::vector<unsigned long> indices_ring1,indices_ring2;
Real d1,d2;
GetIndicesAndMasses(view_ring1, indices_ring1,masses_ring1);
GetIndicesAndMasses(view_ring2, indices_ring2,masses_ring2);
for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
CoordFramePtr frame=traj.GetFrame(i);
d1=frame->GetMinDistBetwCenterOfMassAndView(indices_ring1, masses_ring1, indices_ring2);
d2=frame->GetMinDistBetwCenterOfMassAndView(indices_ring2, masses_ring2, indices_ring1);
dist.push_back(std::min(d1,d2));
}
return dist;
}
void AnalyzeAlphaHelixAxis(const CoordGroupHandle& traj, const EntityView& prot_seg, geom::Vec3List& directions,
geom::Vec3List& centers, unsigned int stride)
//This function calculates the best fitting cylinder to the C-alpha atoms of an EntityView and returns
//the geometric center as well as the axis of that cylinder. We take care to have the axis point towards
//the last residue of the selection, usually the direction of the alpha-helix
{
CheckHandleValidity(traj);
if (!prot_seg.HasAtoms()){
throw Error("EntityView is empty");
}
std::vector<unsigned long> indices_ca;
std::pair<geom::Line3, Real> cylinder;
geom::Line3 axis;
Real sign;
directions.reserve(ceil(traj.GetFrameCount()/float(stride)));
centers.reserve(ceil(traj.GetFrameCount()/float(stride)));
GetCaIndices(prot_seg, indices_ca);
unsigned int n_atoms=indices_ca.size();
for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
CoordFramePtr frame=traj.GetFrame(i);
cylinder=frame->FitCylinder(indices_ca);
axis=cylinder.first;
sign=geom::Dot(axis.GetDirection(),(*frame)[indices_ca[n_atoms-1]]-axis.GetOrigin());
sign=sign/fabs(sign);
directions.push_back(sign*axis.GetDirection());
centers.push_back(axis.GetOrigin());
}
return;
}
//std::vector<geom::Line3> AnalyzeBestFitLine(const CoordGroupHandle& traj, const EntityView& prot_seg,
// unsigned int stride)
void AnalyzeBestFitLine(const CoordGroupHandle& traj, const EntityView& prot_seg, geom::Vec3List& directions,
geom::Vec3List& centers, unsigned int stride)
{
CheckHandleValidity(traj);
if (!prot_seg.HasAtoms()){
throw Error("EntityView is empty");
}
std::vector<unsigned long> indices_ca;
geom::Line3 axis;
directions.reserve(ceil(traj.GetFrameCount()/float(stride)));
centers.reserve(ceil(traj.GetFrameCount()/float(stride)));
GetIndices(prot_seg, indices_ca);
for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
CoordFramePtr frame=traj.GetFrame(i);
axis=frame->GetODRLine(indices_ca);
directions.push_back(axis.GetDirection());
centers.push_back(axis.GetOrigin());
}
return;
}
void AnalyzeBestFitPlane(const CoordGroupHandle& traj, const EntityView& prot_seg, geom::Vec3List& normals,
geom::Vec3List& origins, unsigned int stride)
{
CheckHandleValidity(traj);
if (!prot_seg.HasAtoms()){
throw Error("EntityView is empty");
}
std::vector<unsigned long> indices_ca;
geom::Plane best_plane;
normals.reserve(ceil(traj.GetFrameCount()/float(stride)));
origins.reserve(ceil(traj.GetFrameCount()/float(stride)));
GetIndices(prot_seg, indices_ca);
for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
CoordFramePtr frame=traj.GetFrame(i);
best_plane=frame->GetODRPlane(indices_ca);
normals.push_back(best_plane.GetNormal());
origins.push_back(best_plane.GetOrigin());
}
return;
}
std::vector<Real> AnalyzeHelicity(const CoordGroupHandle& traj, const EntityView& prot_seg,
unsigned int stride)
{
CheckHandleValidity(traj);
if (!prot_seg.HasAtoms()){
throw 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;
}
//This function constructs mean structures from a trajectory
EntityHandle CreateMeanStructure(const CoordGroupHandle& traj, const EntityView& selection,
int from, int to, unsigned int stride)
{
CheckHandleValidity(traj);
if (to==-1)to=traj.GetFrameCount();
unsigned int n_atoms=selection.GetAtomCount();
if (to<from) {
throw Error("to smaller than from");
}
unsigned int n_frames=ceil((to-from)/static_cast<Real>(stride));
if (n_atoms==0){
throw Error("EntityView is empty");
}
if (n_frames<=1) {
throw Error("number of frames is too small");
}
std::vector<unsigned long> indices;
std::vector<geom::Vec3> mean_positions;
EntityHandle eh;
eh=CreateEntityFromView(selection,1);
GetIndices(selection,indices);
mean_positions.assign(n_atoms,geom::Vec3(0.0,0.0,0.0));
for (int i=from; i<to; i+=stride) {
CoordFramePtr frame=traj.GetFrame(i);
for (unsigned int j=0; j<n_atoms; ++j) {
mean_positions[j]+=(*frame)[indices[j]];
}
}
mol::XCSEditor edi=eh.EditXCS(mol::BUFFERED_EDIT);
AtomHandleList atoms=eh.GetAtomList();
for (unsigned int j=0; j<n_atoms; ++j) {
edi.SetAtomPos(atoms[j],mean_positions[j]/Real(n_frames));
}
return eh;
}
Real AnalyzeRMSF(const CoordGroupHandle& traj, const EntityView& selection, int from, int to, unsigned int stride)
// This function extracts the rmsf between two entity views and assigns it
// The views don't have to be from the same entity
// If you want to compare to frame i of the trajectory t, first use t.CopyFrame(i) for example:
// eh=io.LoadPDB(...),t=io.LoadCHARMMTraj(eh,...);Sele=eh.Select(...);t.CopyFrame(0);mol.alg.AnalyzeRMSD(t,Sele,Sele)
{
CheckHandleValidity(traj);
if (to==-1)to=traj.GetFrameCount();
unsigned int n_atoms=selection.GetAtomCount();
if (to<from) {
throw Error("to smaller than from");
}
unsigned int n_frames=ceil((to-from)/static_cast<Real>(stride));
if (n_atoms==0){
throw Error("EntityView is empty");
}
if (n_frames<=1) {
throw Error("number of frames is too small");
}
Real rmsf=0.0;
geom::Vec3 v;
std::vector<unsigned long> sele_indices;
std::vector<geom::Vec3> ref_pos(n_atoms,geom::Vec3(0.,0.,0.));
GetIndices(selection, sele_indices);
for (unsigned int j=0; j<n_atoms; ++j) {
for (int i=from; i<to; i+=stride) {
CoordFramePtr frame=traj.GetFrame(i);
ref_pos[j]+=frame->GetAtomPos(sele_indices[j]);
}
ref_pos[j]/=n_frames;
}
for (int i=from; i<to; i+=stride) {
CoordFramePtr frame=traj.GetFrame(i);
for (unsigned int j=0; j<n_atoms; ++j) {
v=frame->GetAtomPos(sele_indices[j])-ref_pos[j];
rmsf+=geom::Dot(v,v);
}
}
return sqrt(rmsf/float(n_atoms*n_frames));
}
}}} //ns