Commits

Adrien Saladin committed 37e1f21 Merge

merge of cython and heligeom branches

Comments (0)

Files changed (33)

+#include <rigidbody.h>
+#include <pdbio.h>
+#include <iostream>
+#include <BasePair.h>
+#include <Movement.h>
+
+#include "atomselection.h"
+
+using namespace std;
+using namespace PTools;
+
+
+BasePair::BasePair(std::string filename)
+{
+  ReadPDB(filename,rigbody);
+  type = rigbody.GetAtomProperty(0).GetResidType();
+}
+
+
+BasePair::BasePair(const Rigidbody& rigbody)
+{
+  this->rigbody=rigbody;
+  type = rigbody.GetAtomProperty(0).GetResidType();
+}
+
+
+BasePair::~BasePair()
+{
+  
+}
+
+
+string BasePair::printPDB()const
+{
+  return rigbody.PrintPDB ();
+}
+
+std::string BasePair::printPDBofBase(std::string chain) const
+{
+    return rigbody.SelectChainId(chain).CreateRigid().PrintPDB();
+}
+
+void BasePair::setChainID(){
+  unsigned int rigSize=rigbody.Size();
+  for(unsigned int i =0; i< rigSize ; i++)
+  {
+    Atomproperty ap=rigbody.GetAtomProperty(i);
+    if (ap.GetResidType()==type)
+    {
+        ap.SetChainId("A");
+    }
+    else
+    {
+        ap.SetChainId("B");
+    }
+    rigbody.SetAtomProperty(i,ap);
+  }
+}
+
+void BasePair::apply( const Movement& m)
+{
+  m.apply(rigbody);
+}
+
+
+void BasePair::apply(const Matrix& m)
+{
+  apply(Movement (m));
+}
+
+unsigned int BasePair::size() const{
+    return rigbody.Size();
+}
+
+Matrix BasePair::getMatrix() const
+{
+  return rigbody.GetMatrix();
+}
+
+
+Movement BasePair::getMovement()const
+{
+  return Movement(getMatrix());
+}
+
+
+Rigidbody BasePair::getRigidBody()const
+{
+  return rigbody;
+}
+
+
+Rigidbody BasePair::getRigidBodyOfBase(std::string chain)const
+{
+  return rigbody.SelectChainId(chain).CreateRigid();
+}
+
+
+void BasePair::setResID(int idA,int idB)
+{
+  unsigned int baseSize=rigbody.Size();
+  for(unsigned int i =0; i< baseSize ; i++)
+  {
+    Atomproperty ap=rigbody.GetAtomProperty(i);
+    if (ap.GetChainId()=="A")
+    {
+        ap.SetResidId(idA);
+    }
+    else
+    {
+        ap.SetResidId(idB);
+    }     
+    rigbody.SetAtomProperty(i,ap);
+  }
+}
+
+uint BasePair::setAtomNumberOfBase(std::string chain,int num)
+{
+  unsigned int baseSize=rigbody.Size();
+  for(unsigned int i =0; i< baseSize ; i++)
+  {
+    Atomproperty ap=rigbody.GetAtomProperty(i);
+    if (ap.GetChainId()==chain)
+    {
+        ap.SetAtomId(num);
+        num++;
+        rigbody.SetAtomProperty(i,ap);
+    }
+  }
+  return num;
+}
+
+uint BasePair::getResIDofBase(std::string chain)const
+{
+  Atomproperty ap = rigbody.SelectChainId(chain).CreateRigid().GetAtomProperty(0);
+  return ap.GetResidId();
+}
+
+
+void  BasePair::setRigidBody(const Rigidbody& rigbody)
+{
+  this->rigbody=rigbody;
+}
+
+string BasePair::getType() const {
+    return type;
+}
+
+void BasePair::setType(string type) {
+    this->type = type;
+}
+//end namespace
+
+#ifndef BASE_PAIR_H
+#define BASE_PAIR_H
+
+#include <vector>  
+#include <Movement.h>
+
+#include "rigidbody.h"
+
+/*! \brief Rigidbody use by the DNA object
+*
+*/
+namespace PTools
+{
+  
+  class BasePair 
+  {
+    public:
+    ///initialize a new object with a file
+    BasePair(std::string filename);
+    ///initialize a new object from a regular Rigidbody
+    BasePair(const Rigidbody& rigbody);
+    ~BasePair();
+    
+    ///return a string containing the atoms data following the PDB format 
+    std::string printPDB() const;
+
+    ///return a string containing the atoms data of a base (identified by its chain) following the PDB format
+    std::string printPDBofBase(std::string chain) const;
+
+    ///change the chainID of the internal bases to 'A' for the first base (coresponding to the type) and 'B' for the second one
+    void setChainID();
+
+    /// apply a Movement to the BasePair 
+    void apply(const Movement& );
+    /// apply a Matrix to the BasePair
+    void apply(const Matrix&);
+    
+    /// return the Matrix of the BasePair
+    Matrix getMatrix()const;
+    
+    /// return the Movemeny of the BasePair
+    Movement getMovement()const;
+    
+    /// return the Rigidbody of the BasePair
+    Rigidbody getRigidBody()const;
+    /// return the Rigidbody of the specified base (by chain)
+    Rigidbody getRigidBodyOfBase(std::string chain)const;
+    /// define the Rigidbody of the BasePair
+    void setRigidBody(const Rigidbody&);
+    
+    /// return the Residue ID of the specified base (by chain)
+    uint getResIDofBase(std::string chain)const;
+    /// change the ID of res to idA for base on strand A and idB for the one on strand B
+    void setResID(int idA,int idB);
+    /// change the atoms numbers of the specified base (by chain) starting at startnum and returning the last atom number
+    uint setAtomNumberOfBase(std::string chain,int);
+
+    
+    std::string getType() const;
+
+    void setType(std::string type);
+    
+    
+    ///return the number of Atom
+    unsigned int size() const;
+
+    /// return the i-th Atom of the basePair
+    Atom operator[] (uint i) const {
+          if (i>=this->size()) throw std::range_error("DNA: array out of bounds");
+          return rigbody.CopyAtom(i);};
+
+    private:
+    //atribut
+    Rigidbody rigbody;
+    std::string type;
+    
+  };
+  
+}
+
+
+
+#endif // BASE_PAIR_H 
+#include <rigidbody.h>
+#include <pdbio.h>
+#include <superpose.h>
+#include <atomselection.h>
+#include <iostream>
+#include <DNA.h>
+#include <Movement.h>
+#include <fstream>
+#include <limits>
+#include <limits.h>
+#include <numeric>
+#include "Parameter.h"
+
+using namespace std;
+using namespace PTools;
+
+
+//constructor/destructor
+DNA::DNA(string dataBaseFile, string seq, const Movement& mov)
+{
+  //if we want to build the dna from a model
+  if(isPdbFile (seq))
+  {
+      buildDNAfromPDB ( dataBaseFile, seq );
+  }
+  else
+  {
+      assembleSeq (dataBaseFile,seq);
+      applyInitialMov(mov);
+  }
+}
+
+DNA::DNA( const DNA& model )
+{
+    unsigned int modelSize  = model.size();
+    for (uint i =0; i< modelSize; i++)
+    {
+        strand.push_back(model[i]);
+    }
+    
+}
+
+DNA::DNA( string dataBaseFile, Rigidbody model)
+{
+    buildDNAfromModel(dataBaseFile,model);
+}
+
+DNA::DNA()
+{
+    
+}
+void DNA::buildDNAfromPDB (string dataBaseFile, string pdbFile )
+{
+    
+    Rigidbody model = Rigidbody(pdbFile);
+    buildDNAfromModel(dataBaseFile,model);
+}
+
+void DNA::buildDNAfromModel(string dataBaseFile,Rigidbody model)
+{
+    renumberModel (model);
+    model = delSingleBase (model);
+    string seq =getSeq(model);
+    assembleSeq (dataBaseFile,seq);
+
+    placeBasePairs(model);
+}
+
+DNA::~DNA()
+{
+
+}
+
+void DNA::placeBasePairs(const Rigidbody& model)
+{
+    unsigned int DNASize  = (strand.size()*2)-1;
+    unsigned int strandSize  = strand.size();
+    Parameter param =Parameter();
+    for ( unsigned int i = 0; i < strandSize; i++ )// strandSize
+    {
+        Rigidbody modelOfBasePair = getModelOfBasePair( model, i, DNASize-i);
+        strand[i].apply(getMatBetwenBasePair ( modelOfBasePair,i ));
+    }
+}
+
+
+Matrix DNA::getMatBetwenBasePair( const Rigidbody& modelOfBasePair,int pos)const{
+    Parameter param =Parameter();
+    Superpose_t sup =superpose (param.buildAxisCentered(modelOfBasePair),param.buildAxisCentered(strand[pos].getRigidBody()));
+    return sup.matrix;
+}
+
+
+void DNA::renumberModel (Rigidbody& model)const
+{
+
+  unsigned int tempId=  model.GetAtomProperty(0).GetResidId();
+  string chain;
+  unsigned int nbRes=0;
+  unsigned int second = 0;
+
+    unsigned int strandSize = 0;
+    if ((model.SelectAtomType("C1'").Size() / 2) > 0) {
+        strandSize = model.SelectAtomType("C1'").Size() / 2;
+    } else if ((model.SelectAtomType("C1*").Size() / 2) > 0) {
+        strandSize = model.SelectAtomType("C1*").Size() / 2;
+    } else if((model.SelectAtomType("GS1").Size() / 2) > 0) {
+        strandSize = model.SelectAtomType("GS1").Size() / 2;
+    }
+
+  bool isjumna = isJumna(model);
+  unsigned int modelSize=model.Size();
+  chain = "A";
+  for (unsigned int i =0; i < modelSize; i++ )
+  {
+    Atomproperty ap=model.GetAtomProperty(i);
+    unsigned int Id = ap.GetResidId();
+    if ( tempId != Id )
+    {
+        
+        if (nbRes >= strandSize -1){
+            chain = "B";
+            if (isjumna)
+            {
+                if (second == 0)
+                {
+                    second=nbRes*2 +1;
+                }
+                nbRes= second--;
+                tempId =Id;
+            }
+            else
+            {
+                tempId =Id;
+                nbRes++;
+            }
+        }
+        else
+        {
+            tempId =Id;
+            nbRes++;
+        }
+    }
+    ap.SetResidId(nbRes);
+    ap.SetChainId(chain);
+    model.SetAtomProperty(i,ap);
+  }
+}
+
+bool DNA::isJumna (const Rigidbody& model)const
+{
+   //jumna have an inversed numerotation so on a 10 base dna (strand A: 0,1,2,3,4 and strand B:5,6,7,8,9)
+   //the base 0 is associated with 5 instead of 9.
+   //to check wich convention is followed, I'm gonna test the distance between 0-5 and 0-9, the shorter being the correct coupling
+    AtomSelection sel = model.SelectAtomType("C1'");
+    if (sel.Size() == 0) sel = model.SelectAtomType("C1*");
+    if (sel.Size() == 0) sel = model.SelectAtomType("GS1");
+
+
+    double d1 = Dist(sel[0],sel[sel.Size()-1]);
+    double d2 = Dist(sel[0],sel[sel.Size()/2]);
+
+    return (d1>d2);
+}
+
+Rigidbody DNA::delSingleBase (Rigidbody& model)const
+{
+    string seq;
+    unsigned int strandSize;
+    Rigidbody newModel = Rigidbody();
+    
+    if ((model.SelectAtomType("C1'").Size()) >0)
+    {
+        strandSize = model.SelectAtomType("C1'").Size();
+    }else if ((model.SelectAtomType("C1*").Size()) >0)
+    {
+        strandSize = model.SelectAtomType("C1*").Size();
+    }else if ((model.SelectAtomType("GS1").Size()) >0)
+    {
+        strandSize = model.SelectAtomType("GS1").Size();
+    }else {return model;}
+    for ( unsigned int i=0 ; i< strandSize ; i++ )
+    {
+         string type = model.SelectResRange( i, i)[0].GetResidType();
+         // /!\ the order of the check is important! somme pdb use a CYT description for C, a wrong order could detect this as a T
+         if      ( type.find ('G') != string::npos || type.find ('g') != string::npos) seq+='G';
+         else if ( type.find ('C') != string::npos || type.find ('c') != string::npos) seq+='C';
+         else if ( type.find ('A') != string::npos || type.find ('a') != string::npos) seq+='A';
+         else if ( type.find ('T') != string::npos || type.find ('t') != string::npos) seq+='T';
+         else {
+            cerr <<"unrecognised resid type for base "<< i+1 <<" : "<<type<< endl;
+            exit(0);
+         };
+
+    }
+    int l = seq.length();
+    
+    int solution[l];
+    for (int i=0; i<l;i++){solution[i]=0;}
+
+    string rseq = seq;
+    for (int i = 0; i<l;i++){
+        if      (seq[l-i-1]=='A')  rseq[i]='T';
+        else if (seq[l-i-1]=='T')  rseq[i]='A';
+        else if (seq[l-i-1]=='G')  rseq[i]='C';
+        else if (seq[l-i-1]=='C')  rseq[i]='G';
+    }
+    //cout << seq << " " << rseq.length() << endl;
+    if (isAlign(seq,rseq))
+    {
+        return model;
+    }
+    else
+    {
+        //build matrix
+	int mat [l][l];
+        for (int i =0; i<l; i++){
+            for(int j =0; j<l; j++){
+                mat[i][j] = 0;
+            }
+        }
+	//fill it
+        for (int i=l-1; i>=0;i--)
+        {
+            for(int j =0; j<l; j++)
+            {
+                if ((i == j) || (i-1 == j)|| (i>j))
+                {
+                    mat[i][j] = 0;
+                }
+                else
+                {
+                    int one =0;
+                    if ((seq[i]=='A' && seq[j]=='T') || (seq[i]=='T' && seq[j]=='A') || (seq[i]=='G' && seq[j]=='C') || (seq[i]=='C' && seq[j]=='G'))
+                    {
+                         one = mat[i+1][j-1] +1;
+                    }
+                    int two = mat[i+1][j];
+		    int three=mat[i][j-1];
+		    mat[i][j] = max(max(one, two),three);
+                  
+                }
+            }
+        }
+	//go trougth it to find aligned part
+
+        int i = 0;
+        int j= l-1;
+        int current =mat[i][j];
+        while (current !=0)
+        {
+            int pair = mat[i+1][j-1];
+            if(pair < mat[i+1][j])
+            {
+                current=mat[i+1][j];
+                i++;
+            }
+            else if (pair< mat[i][j-1])
+            {
+                current=mat[i][j-1];
+		j--;
+            }
+            else
+            {
+                if ((seq[i]=='A' && seq[j]=='T') || (seq[i]=='T' && seq[j]=='A') || (seq[i]=='G' && seq[j]=='C') || (seq[i]=='C' && seq[j]=='G'))
+                {
+			solution[i]=1;
+			solution[j]=1;
+                }
+		i++;
+		j--;
+		current = pair;
+            }
+        }
+
+
+        //print matrix
+      /*  for (int i =0; i<l; i++){
+            for(int j =0; j<l; j++){
+                cout << mat[i][j];
+            }
+            cout << endl;
+        }
+    cout <<  endl;
+    cout <<  endl;
+    */
+    //check ambiguity
+    for (int i =0; i<l; i++)
+    {
+        if(solution[i]==0)
+        {
+            char c = seq[i];
+            int start = i;
+            int end = i;
+            //find the start
+            while (seq[start]==c) start--;
+            start++;
+            while (seq[end]==c) end++;
+            end--;
+            if ((start-end) == 0) continue;
+            int diff [end-start];
+            for (int j =start;j<end;j++)
+            {
+                diff[j-start]= calcPart (solution,j,i,l ); 
+            }
+            int min = diff[0];
+            int idx = start;
+            for (int j =1; j<end; j ++)
+            {
+                if (diff[j]<min)
+                {
+                    min = diff[j];
+                    idx = j+start;
+                }
+            }
+            solution[i] = 1;
+            solution[idx] = 0;
+            
+        }
+    }
+    //for (int i =0; i<l; i++)
+    //{
+    //    cout << solution[i];
+    //}
+    //cout << endl;
+ 
+
+  
+    for (int i =0; i<l; i++)
+    {
+        if ( solution[i]==1)
+        {
+              newModel= newModel + model.SelectResRange(i,i).CreateRigid();
+        }
+    }
+    //cout << newModel.PrintPDB()<< endl;
+    renumberModel (newModel);
+    
+
+    return newModel;
+    }
+}
+
+int DNA::calcPart (int solution[],int pos0, int pos1, int l )const
+{
+    int alt[l];
+    for (int j =0; j<l; j++)
+    {
+        alt[j] = solution[j];
+    }
+    alt [pos1] = 1;
+    alt [pos0] = 0;
+
+    
+    int firstAlt = 0;
+    int secondAlt = -1;
+
+    for(int j =0; j<l; j++)
+    {
+       if (alt[j] == 1 && secondAlt == -1)
+       {
+           firstAlt ++;
+       }
+       if (alt[j] ==0 && firstAlt != 0)
+       {
+           secondAlt = 0;
+       }
+       if (alt[j] == 1 && secondAlt != -1)
+       {
+           secondAlt +=1;
+       }
+    }
+    return abs(firstAlt-secondAlt);
+}
+
+bool DNA::isAlign(std::string s1,std::string s2,int shift)const
+{
+    return (s1.compare(shift,s1.length(), s2, 0,s2.length()-shift)==0);
+}
+
+Rigidbody DNA::getModelOfBasePair(const Rigidbody& model,int posA,int posB)const
+{
+        return (model.SelectResRange(posA, posA)|model.SelectResRange(posB, posB)).CreateRigid();
+
+
+}
+
+
+void DNA::assembleSeq (std::string dataBaseFile, std::string seq)
+{
+      Rigidbody dataBase = Rigidbody(dataBaseFile);      
+      string chainIDs = getChainIDs(dataBase);
+
+      //"map" for the rigidbody, an iD corespond to its rigidbody
+      vector<Rigidbody> vbase = buildVbase(chainIDs,dataBase);
+      //build the strand from the seq
+      buildStrand(seq, chainIDs, vbase);
+      //make sure that every BasePaire have a different id
+      changeFormat();
+}
+
+string DNA::getSeq ( const Rigidbody& model)const
+{
+    string seq;
+    unsigned int strandSize;
+    if ((model.SelectAtomType("C1'").Size()/2) >0)
+    {
+        strandSize = model.SelectAtomType("C1'").Size()/2;
+    }else if ((model.SelectAtomType("C1*").Size()/2) >0)
+    {
+        strandSize = model.SelectAtomType("C1*").Size()/2;
+    }else if ((model.SelectAtomType("GS1").Size()/2) >0)
+    {
+        strandSize = model.SelectAtomType("GS1").Size()/2;
+    }else {return "";}
+    for ( unsigned int i=0 ; i< strandSize ; i++ )
+    {
+         string type = model.SelectResRange( i, i)[0].GetResidType();
+         // /!\ the order of the check is important! somme pdb use a CYT description for C, a wrong order could detect this as a T
+         if      ( type.find ('G') != string::npos || type.find ('g') != string::npos) seq+='G';
+         else if ( type.find ('C') != string::npos || type.find ('c') != string::npos) seq+='C';
+         else if ( type.find ('A') != string::npos || type.find ('a') != string::npos) seq+='A';
+         else if ( type.find ('T') != string::npos || type.find ('t') != string::npos) seq+='T';
+
+    }
+    return seq;
+}
+
+
+vector<Rigidbody> DNA::buildVbase(string chainIDs, Rigidbody& dataBase)const
+{
+  vector<Rigidbody> vbase;
+  unsigned int chainIDsSize = chainIDs.size();
+  for (unsigned int i = 0; i < chainIDsSize ; i++)
+  {
+
+    vbase.push_back(dataBase.SelectChainId(chainIDs.substr(i,1)).CreateRigid());
+  }
+  return vbase;
+}
+
+
+string DNA::getChainIDs(const Rigidbody& rb)const
+{
+  string chainIDs;
+  AtomSelection selection = rb.SelectAllAtoms ();
+  string tmp = "";
+  unsigned int selectionSize = selection.Size();
+
+  for (unsigned int i=0; i < selectionSize ;i++)
+  {
+    string id = selection[i].GetChainId();
+    if(id !=tmp)
+    {
+      tmp=id;
+      chainIDs += id;
+    }
+  }
+  return chainIDs;
+}
+
+bool DNA::isPdbFile (std::string seq) const
+{
+    return ( (seq.size() >=3) && ((seq.substr(seq.size()-3,seq.size())=="pdb") ||(seq.substr(seq.size()-3,seq.size())=="red")) );
+}
+
+void DNA::buildStrand(std::string seq, std::string chainIDs, const std::vector<Rigidbody>& vbase)
+{
+  unsigned int seqSize = seq.size();
+  unsigned int chainIDsSize = chainIDs.size();
+
+  for (unsigned int i =0; i < seqSize; i++ )
+  {
+    for (unsigned int j =0; j < chainIDsSize; j++ )
+    {     
+      if (seq[i] == chainIDs[j])
+      {
+	strand.push_back(BasePair(vbase[j]));
+      }
+    }    
+  }  
+}
+
+
+void DNA::changeFormat()
+{
+  unsigned int nbAtom = 0;
+  unsigned int strandSize  = strand.size();
+  for (unsigned int i =0; i < strandSize; i++ )
+  {
+    //corect ID chain
+    strand[i].setChainID();
+    //numerotation residu
+    strand[i].setResID(i,(strandSize + strandSize-1)-i);
+    //numerotation atom (first part)
+    nbAtom = strand[i].setAtomNumberOfBase("A",nbAtom);
+
+  }
+  for (unsigned int i = strandSize-1; i > 0 ; i-- )
+  {
+    //numerotation atom (second part)
+    nbAtom = strand[i].setAtomNumberOfBase("B",nbAtom);
+  }
+  //last part of atom numerotation (because of error with decreasing unsigned int)
+  strand[0].setAtomNumberOfBase("B",nbAtom);
+}
+
+
+void DNA::applyInitialMov(const Movement& mov)
+{
+  unsigned int strandSize  = strand.size();
+  for (unsigned int i=1; i <strandSize; i++)
+  {
+    strand[i].apply(strand[i-1].getMovement());
+    strand[i].apply(mov);
+  }
+}
+
+
+unsigned int DNA::size()const
+{
+  return strand.size();
+}
+
+unsigned int DNA::atomSize() const
+{
+    uint tot=0;
+    uint strandSize= this->size();
+    for (uint i=0; i<strandSize;i++)
+    {
+        tot+=strand[i].size();
+    }
+    return tot;
+}
+
+
+string DNA::printPDB()const
+{
+  string strandA, strandB;
+  unsigned int strandSize  = strand.size();
+  for ( unsigned int i =0; i < strandSize ; i++ )
+  {
+    strandA += strand[i].printPDBofBase("A");
+    strandB += strand[strandSize-1-i].printPDBofBase("B");
+  }
+  string out= strandA + strandB;
+  return out.substr(0,out.size()-1);
+}
+
+std::string DNA::printPDBofStrand( std::string chain ) const
+{
+  string out;
+  unsigned int strandSize  = strand.size();
+  for ( unsigned int i =0; i < strandSize ; i++ )
+  {
+    out += strand[i].printPDBofBase( chain );
+  }
+  return out.substr(0,out.size()-1);
+}
+
+string DNA::printParam() const
+{
+  stringstream ss;
+  unsigned int strandSize  = strand.size();
+  for ( unsigned int i =1; i < strandSize ; i++ )
+  {
+    ss << "BasePair "<< i-1 <<"->"<<i<< " : "<<getLocalParameter(i).toFormatedString()+"\n";
+  }
+  return ss.str().substr(0,ss.str().size()-1);
+}
+
+
+void DNA::writePDB(std::string fileName)const
+{
+  ofstream myfile;
+  myfile.open(fileName.c_str());
+  myfile << printPDB();
+  myfile.close();
+}
+
+Rigidbody DNA::createRigid()const
+{
+  //we use the method createRigidOfStrand() to make sure that the rigidbody returned follow the established format
+  Rigidbody chainA = createRigidOfStrand("A");
+  Rigidbody chainB = createRigidOfStrand("B");
+
+  //we add the atoms of the chain B to the chainA
+  unsigned int chainBsize  = chainB.Size();
+  for ( unsigned int i =0; i < chainBsize ; i++ )
+  {
+      chainA.AddAtom(chainB.CopyAtom(i));
+  }
+  return chainA;
+}
+
+Rigidbody DNA::createRigidOfStrand(std::string chain)const
+{
+      Rigidbody dna;
+  unsigned int strandSize  = strand.size();
+
+
+  for ( unsigned int i =0; i < strandSize ; i++ )
+  {
+      Rigidbody basePair;
+      if (chain == "B" || chain == "b")
+      {
+          basePair = strand[strandSize-1-i].getRigidBodyOfBase(chain);
+      }
+      else
+      {
+          basePair = strand[i].getRigidBodyOfBase(chain);
+      }
+      unsigned int basePairSize  = basePair.Size();
+      for (unsigned int j=0; j <basePairSize ; j++)
+      {
+          Atom a =basePair.CopyAtom(j);
+          dna.AddAtom(a);
+      }
+  }
+  return dna;
+}
+
+void DNA::changeRepresentation(std::string dataBaseFile)
+{
+  Rigidbody dataBase = Rigidbody(dataBaseFile);
+  string chainIDs = getChainIDs(dataBase);
+  
+  //"map" for the rigidbody, an iD corespond to its rigidbody
+  vector<Rigidbody> vbase = buildVbase(chainIDs,dataBase);
+
+    unsigned int chainIDsSize = chainIDs.size();
+
+  unsigned int strandSize  = strand.size();
+  for (unsigned int i = 0; i < strandSize ; i++)
+  {
+    Movement mov = Movement(strand[i].getMatrix());
+
+    for (unsigned int j =0; j < chainIDsSize; j++ )
+    {     
+      if ( strand[i].getType()[0] == chainIDs[j])
+      {
+	strand[i]=BasePair(vbase[j]);
+	strand[i].apply(mov);
+      }
+    }
+  }
+  
+  changeFormat();
+}
+
+
+Matrix DNA::getLocalMatrix(int pos)const
+{
+   if (pos == 0) return strand[0].getMatrix();
+  
+   Matrix mtot = strand[pos].getMatrix();
+   Matrix prec = inverseMatrix44(strand[pos-1].getMatrix());
+
+   return matrixMultiply( prec, mtot );
+   
+}
+
+Parameter DNA::getLocalParameter(int pos)const
+{
+  Parameter param =Parameter();
+  param.MeasureParameters(param.buildAxisCentered(strand[pos-1].getRigidBody()),param.buildAxisCentered(strand[pos].getRigidBody()));
+  return param;
+
+}
+
+void DNA::applylocalMov(const Movement& mov,int pos)
+{
+  Matrix nextlocal = getLocalMatrix(pos+1);
+  strand[pos].apply(mov);
+
+  unsigned int strandSize  = strand.size();
+  for (unsigned int i=pos+1; i <strandSize; i++)
+  {
+    nextlocal = reconstruct(i,nextlocal);
+  }
+}
+
+
+void DNA::applyglobalMov(const Movement& mov)
+{
+  Matrix nextlocal;
+  unsigned int strandSize  = strand.size();
+  if (strandSize>1){
+    nextlocal = getLocalMatrix(1);
+  }
+  strand[0].apply(mov);
+
+  for (unsigned int i=1; i <strandSize; i++)
+  {
+    nextlocal = reconstruct(i,nextlocal);
+    strand[i].apply(mov);
+  }  
+    
+}
+
+
+void DNA::applyLocal(const Movement& mov,int posMov, int posAnchor)
+{
+  BasePair anchor = strand[posAnchor];
+  applylocalMov(mov,posMov);
+  relocate(anchor,posAnchor);
+}
+
+
+void DNA::applyLocal(const Matrix&  m,int posMov, int posAnchor)
+{
+  applyLocal ( Movement(m), posMov, posAnchor);
+}
+
+
+void DNA::applyGlobal(const Movement& mov ,int posAnchor)
+{
+  BasePair anchor = strand[posAnchor];
+  applyglobalMov(mov);
+  relocate(anchor,posAnchor);
+}
+
+
+void DNA::applyGlobal(const Matrix& m ,int posAnchor)
+{
+  applyLocal ( Movement(m), posAnchor);
+}
+
+
+void DNA::apply(const Movement& mov)
+{
+    DNA::apply(mov.getMatrix());
+}
+
+
+void DNA::apply(const Matrix& m)
+{
+  unsigned int strandSize  = strand.size();
+  for (unsigned int i=0; i <strandSize; i++)
+  {
+    Rigidbody rb = strand[i].getRigidBody();
+    rb.ApplyMatrix(m);
+    strand[i].setRigidBody(rb);
+  }
+}
+
+
+double DNA::Rmsd(const DNA& model)const
+{
+
+    double aSize=atomSize();
+    if (model.atomSize() != aSize)
+    {
+        std::cerr << "Error: trying to superpose two DNA of different sizes" << std::endl ;
+        throw std::invalid_argument("RmsdSizesDiffers");
+    }
+
+    double total = 0.0;
+    uint strandSize = size();
+    for (uint i = 0; i<strandSize;i++){
+        uint pbSize= strand[i].size();
+        for (uint j =0; j<pbSize;j++)
+        {
+                Atom atom1=model[i][j];
+		Atom atom2=strand[i][j];
+
+		total+=Dist2(atom1,atom2);
+        }
+    }
+    return sqrt(total/(dbl) aSize) ;
+
+}
+
+void DNA::changeBasePair(const BasePair& bp, int pos)
+{
+    strand[pos]=bp;
+}
+
+void DNA::relocate(const BasePair& anchor,int posAnchor)
+{
+  apply(superpose(anchor.getRigidBody(),strand[posAnchor].getRigidBody(),0).matrix);
+}
+
+
+Matrix DNA::reconstruct(int pos, const Matrix& local)
+{
+  Matrix nextlocal;
+  if ((unsigned int)pos<strand.size()) nextlocal = getLocalMatrix(pos+1);
+  
+  Matrix prec = strand[pos-1].getMatrix();
+  
+  strand[pos].apply(inverseMatrix44(strand[pos].getMatrix()));//erasing the old matrix
+  
+  strand[pos].apply(prec);
+  strand[pos].apply(local);
+  
+  return nextlocal;
+}
+
+
+void DNA::add(const DNA & d, const Movement & mov)
+{
+    this->add(d[0], mov);
+    for(uint i =1; i< d.size();i++)
+    {
+        this->add(d[i],Movement(d.getLocalMatrix(i)));
+    }
+}
+
+/// add a basepair at the end of the strand of this DNA
+void DNA::add(BasePair bp, const Movement & mov)
+{
+    Matrix oldmouvement = bp.getMatrix ();
+    bp.apply(inverseMatrix44 (oldmouvement));
+    if (strand.size()>0)
+    {
+        bp.apply(matrixMultiply(strand[strand.size()-1].getMatrix(),mov.getMatrix()));
+    }
+    else bp.apply(mov.getMatrix());
+    strand.push_back(bp);
+    this->changeFormat();
+}
+
+DNA DNA::subDNA(int start,int end)const
+{
+    DNA newdna = DNA();
+
+    newdna.add(strand[start],Movement(strand[start].getMatrix()));
+    for (int i=start+1;i<end;i++)
+    {
+        newdna.add(strand[i],Movement(this->getLocalMatrix(i)));
+    }
+    return newdna;
+}
+
+void DNA::replace(const DNA & d,int start)
+{
+    DNA preDNA = this->subDNA(0,start);
+    DNA postDNA =this->subDNA(start+d.size(),this->size());
+
+    Movement initMov = Movement(strand[0].getMatrix());
+    strand.clear();
+
+
+    this->add(preDNA,initMov);
+    this->add(d);
+    this->add(postDNA);
+
+
+}
+
+void DNA::changeType(int pos, std::string type, std::string filename) {
+    Rigidbody dataBase = Rigidbody(filename);
+    Movement mov = Movement(strand[pos].getMatrix());
+
+    strand[pos] = BasePair(dataBase.SelectChainId(type).CreateRigid());
+    strand[pos].apply(mov);
+
+    changeFormat();
+}
+
+void DNA::translate(Coord3D coord)
+{
+  unsigned int strandSize  = strand.size();
+  for (unsigned int i=0; i <strandSize; i++)
+  {
+    Rigidbody rb = strand[i].getRigidBody();
+    rb.Translate(coord);
+    strand[i].setRigidBody(rb);
+  }
+
+}
+
+#ifndef DNA_H
+#define DNA_H
+
+#include <vector>  
+#include <BasePair.h>
+#include <Movement.h>
+#include <Parameter.h>
+#include "rigidbody.h"
+
+namespace PTools
+{
+  
+  class DNA 
+  {
+    
+    public:
+    ///initialize a new object with a sequence and a database of pdb to pick from. an initial movement for the construction of the dna can be precised. 
+    DNA( std::string , std::string , const Movement & mov = BDNA());
+    ///initialize a dna from another dna.
+    DNA( const DNA& model );
+    ///initialize a dna from dna in a rigidbody
+    DNA( std::string dataBaseFile,Rigidbody model);
+    ///initialize a hollow dna (no BasePair in the strand)
+    DNA();
+    ~DNA();
+
+
+    ///return the number of BasePair
+    unsigned int size() const;
+    
+    ///return the number of Atom
+    unsigned int atomSize() const;
+
+    ///return a string containing the atoms data following the PDB format
+    std::string printPDB() const;
+
+    ///return a string containing the atoms data following the PDB format
+    std::string printPDBofStrand( std::string chain ) const;
+
+    ///return a string containing all the local parameter formated
+    std::string printParam() const;
+
+    ///write in a file the atoms data following the PDB format
+    void writePDB(std::string)const;
+    
+    ///change the database use for the BasePair. the new database must contain the same names for pdb that the old one.
+    void changeRepresentation(std::string);
+    
+    ///apply a Movement to a specified BasePair. you can specify an anchor 
+    void applyLocal(const Movement&,int posMov, int posAnchor = 0);
+    ///apply a Movement to all the BasePairs and reposition the DNA according to the anchor
+    void applyGlobal(const Movement&,int posAnchor);
+    ///apply a Movement to the DNA as a rigidbody
+    void apply(const Movement&);
+
+    ///apply a Matrix to a specified BasePair. you can specify an anchor 
+    void applyLocal(const Matrix&,int posMov, int posAnchor = 0);
+    ///apply a Matrix to all the BasePairs and reposition the DNA according to the anchor
+    void applyGlobal(const Matrix&,int posAnchor);
+    ///apply a Matrix to the DNA as a rigidbody
+    void apply(const Matrix&);
+
+    ///apply a vector to the DNA as a rigidbody
+    void translate(Coord3D coord);
+    
+    ///return the local Matrix of the specified BasePair (for the position i the Matrix to go from i-1 to i)
+    Matrix getLocalMatrix(int pos)const;
+    ///return the local Parameter of the specified BasePair (for the position i the Matrix to go from i-1 to i)
+    Parameter getLocalParameter(int pos)const;
+
+    ///return a Rigidbody of the DNA()
+    Rigidbody createRigid()const;
+    ///return a Rigidbody of the strand
+    Rigidbody createRigidOfStrand(std::string chain)const;
+
+    
+    //replace the base pair at the indicated position by the new base pair
+    void changeBasePair(const BasePair& bp, int pos);
+
+
+    /// return the RMSD between two DNA (not aligning them and assuming they are comparable)
+    double Rmsd(const DNA&)const;
+
+    ///------edition functions------
+
+    /// return the i-th BasePair of the strand
+    BasePair operator[] (uint i) const {
+          if (i>=this->size()) throw std::range_error("DNA: array out of bounds");
+          return strand[i];};
+
+    /// add the basePairs of a DNA to the strand of this DNA. the specified movement do the liason betwen the two strand
+    void add(const DNA & d, const Movement & mov = BDNA());
+    /// add a basepair at the end of the strand of this DNA
+    void add(BasePair bp, const Movement & mov = BDNA());
+
+    ///return the specified subDNA
+    DNA subDNA(int start,int end)const;
+
+    ///replace the basePair of this DNA by the basePair of a given DNA starting from a given position to the end of one of the DNAs
+    void replace(const DNA & d,int start);
+
+    ///change the type of the base pair at position pos to the indicated type using the coresponding structure in the designed database of pdb.
+    void changeType(int pos, std::string type, std::string filename );
+
+
+
+
+    private:
+
+    //attribut
+    std::vector<BasePair> strand;
+    
+    //methods
+    ///return a string with the seq of all the diferent chain ID of a Rigidbody.
+    std::string getChainIDs(const Rigidbody&)const;
+    ///construct the strand from the sequence, the  chain Ids and a corresponding vector of Rigidbody
+    void buildStrand(std::string,std::string,const std::vector<Rigidbody> &);
+    /// apply an initial Movement during the initialisation of the DNA
+    void applyInitialMov(const Movement&);
+    ///rebuild the DNA from a specified position. useful for taking into account a change in a base.
+    Matrix reconstruct(int pos,const Matrix&);
+    ///give the BasePair a correct number according to they rank in the strand
+    void changeFormat();
+    ///apply a Movement to a specified BasePair.
+    void applylocalMov(const Movement&,int pos);
+    ///apply a Movement to all the BasePairs
+    void applyglobalMov(const Movement& );
+    ///reposition the DNA according to the anchor
+    void relocate(const BasePair& anchor,int posAnchor);
+
+    ///initialize a new object from a PDB file
+    void buildDNAfromPDB ( std::string database, std::string pdbfile);
+    void buildDNAfromModel(std::string dataBaseFile,Rigidbody model  );
+    ///check if a file have a pdb extension
+    bool isPdbFile (std::string)const;
+    ///return the sequence of a DNA model
+    std::string getSeq (const Rigidbody& model)const;
+    ///return a vector with an example of each base
+    std::vector<Rigidbody> buildVbase(std::string chainIDs, Rigidbody& dataBase)const;
+    ///create the different base in strand following the seq order.
+    void assembleSeq (std::string dataBaseFile, std::string seq);
+
+    ///change the numerotation of the base in the model to conform to the following order
+    ///A 1 2 3
+    ///B 6 5 4
+    void renumberModel (Rigidbody& model)const;
+   
+
+    // test if the model follow the jumna convention
+    //A 1 2 3
+    //B 4 5 6
+    bool isJumna (const Rigidbody& model)const;
+
+    //detect and delete single base at extremeties of DNA
+    Rigidbody delSingleBase (Rigidbody& model)const;
+
+    //use by delsinglebase
+    int calcPart (int solution [] , int pos1, int pos2, int l )const;
+
+    //detect if two string are aligned. a shift can be precised
+    bool isAlign(std::string s1,std::string s2,int shift = 0)const;
+    
+    ///place the base in the same position as the model
+    void placeBasePairs(const Rigidbody& model);
+
+    ///return the base pair composed of the base on posA in the strand A of model and the base on posB in the strand B of model
+    Rigidbody getModelOfBasePair(const Rigidbody& model,int posA,int posB)const;
+
+    ///return the matrix betwen a model base pair and the basepair on pos in strand
+    Matrix getMatBetwenBasePair(const Rigidbody& modelOfBasePair,int pos)const;
+
+
+  };
+  
+}//end namespace
+
+#endif // DNA_H 

DNA/buildDNAalongAnAxis.py

+from ptools import *
+import math
+import sys
+import os
+
+PB_CG = "pb.cg.pdb"
+PB_AA = "pb.aa.pdb"
+
+def readline(lineName):
+    line=[]
+    l = Rigidbody(lineName)
+    for i in xrange(l.Size()):
+        line.append(l.CopyAtom(i))
+    return line
+    
+def cleanPDB(rig):
+    rSize = rig.Size()
+    l=[]
+    base =[]
+    nb=-1
+    oldChain = None
+    for i in xrange(0,rSize):
+        a=rig.CopyAtom(i)
+        chain = a.GetChainId ()
+        resid = a.GetResidId ()
+        if nb != resid:
+            if nb != -1:
+                l.append([oldChain,base])
+            nb = resid
+            base=[]
+        if oldChain != chain:
+            oldChain = chain
+        
+        base.append(a)
+    l.append([chain,base])
+    
+    l = sorted(l, key=lambda base: base[0])
+    
+    for i,b in enumerate(l):
+        if b[0] == "B":
+            break
+    l[i],l[i+1]=l[i+1],l[i]
+    
+    newRig = Rigidbody()
+    nbAtom=0
+    for nbRes,base in enumerate(l):
+        #base = base[1]
+        
+        for a in base[1]:
+            a.SetResidId (nbRes)
+            a.SetAtomId (nbAtom)
+            nbAtom += 1
+            newRig.AddAtom(a)
+            #print a.ToPdbString(),
+    return newRig
+
+
+def ScalProd(a,b):
+    return a.x * b.x + a.y * b.y + a.z * b.z ;
+
+def nextModel(line,D,pos):
+    start = line [pos]
+    nextpos = None
+    for i,dot in enumerate(line[pos+1:]):
+        if Dist (dot,start)>D:
+            nextpos=pos+1+i-1
+            A= line[pos+1+i-1]
+            B=line[pos+1+i]
+            break
+    if not nextpos:
+        return None,pos
+    AB = B.GetCoords () - A.GetCoords()
+    Astart = start.GetCoords () - A.GetCoords()
+    AB.Normalize()
+    dproj = ScalProd (AB,Astart)
+    proj = Atom(Atomproperty (),A.GetCoords()+(AB*dproj))
+    startProj = Dist(start,proj)
+    distX = math.sqrt((D**2)-(startProj**2))#V
+    at = Atomproperty ()
+    at.SetAtomId (pos)
+    X= Atom(at,proj.GetCoords()+(AB*distX))
+
+    
+    model =Rigidbody()
+    model.AddAtom ( start )
+    model.AddAtom ( X )
+    return model,nextpos
+        
+def buildFirstPB(mobil,D):
+    dna = DNA(PB_CG ,"A",BDNA())
+    p = Parameter ()
+    axe = p.buildAxisCGGeometricCenter (dna[0].getRigidBody ())
+    axeZ= Rigidbody()
+    axeZ.AddAtom (axe.CopyAtom (0))
+
+    axeZ.AddAtom ( Atom(Atomproperty (),axe.CopyAtom (3).GetCoords()*abs(D)))
+    m = superpose(mobil,axeZ).matrix
+    #axeZ.apply(m)
+    #print axeZ.PrintPDB()
+    #print mobil.PrintPDB()
+    #print dna.printPDB()
+    dna.apply(m)
+    #print dna.printPDB()
+    return dna
+
+def buildNextPB(dna,mobil,D):
+    new = DNA(PB_CG ,"A",BDNA())
+    #print mobil.PrintPDB()
+    p = Parameter ()
+    axe = p.buildAxisCGGeometricCenter (dna[dna.size()-1].getRigidBody ())
+    
+    
+    model = Rigidbody()
+    model.AddAtom (mobil.CopyAtom (1)) 
+    
+    theo = Rigidbody()
+    theo.AddAtom (Atom(Atomproperty(),axe.CopyAtom (0).GetCoords() + (axe.CopyAtom (3).GetCoords() - axe.CopyAtom (0).GetCoords() ).Normalize ()*abs(D)))
+    
+    m = superpose(model,theo).matrix
+    
+    mov = Movement (m)+BDNA()
+    dna.add(new,mov)
+    return dna
+
+def buildDNA(line):
+    
+    dnaRig = Rigidbody()
+
+    segment = DNA(PB_AA ,"AA",BDNA())
+    center0 =Atom(Atomproperty (), segment[0].getRigidBody ().FindCenter ())
+    center1 =Atom(Atomproperty (), segment[1].getRigidBody ().FindCenter ())
+    D = Dist(center0,center1)
+    pos = 0
+    mobil=Rigidbody()
+    mobil.AddAtom ( center0 )
+    mobil.AddAtom ( center1 )
+   
+    
+    model,newpos = nextModel(line,D,pos)
+    m = superpose(model,mobil).matrix
+    
+    
+    segment.apply(m)
+    mobil.ApplyMatrix(m)
+    
+    dna = buildFirstPB(mobil,D)
+    dna = buildNextPB(dna,mobil,D)
+    dna.changeRepresentation(PB_AA)
+
+    segment = DNA(dna)
+    #print dna.createRigid().PrintPDB()
+    dnaRig= dnaRig + segment.createRigid()
+    #dnaRig= dnaRig + dna.createRigid()
+    i=0
+    size = 2
+    while True :
+       
+        pos = newpos
+        
+        model,newpos = nextModel(line,D,pos)
+        if newpos == pos:
+            break
+        m = superpose(model,mobil).matrix
+        mobil.ApplyMatrix(m)
+        segment.apply(m)
+        m=(Twist( 35.9063052632 )+Roll( -2.66592947368 )+Tilt( -1.80234789474 )+Slide( -1.34487389474 )+Shift( -0.425181378947 )).getMatrix();
+        newPB=segment[1]
+        
+        i+=1
+        for i in xrange(0,size-1):
+            newPB.apply(m)
+            
+        #dna = buildNextPB(dna,mobil,D)
+        dnaRig= dnaRig + newPB.getRigidBody ()
+        size+=1
+
+    #print dna.printPDB()
+    dnaRig=cleanPDB(dnaRig)
+    return DNA(PB_AA,dnaRig)
+        
+    
+
+def main():
+    nargs = len(sys.argv)
+    if nargs < 2:
+        print "usage: buildDNAalongAnAxis.py  axis"
+        raise SystemExit
+
+    line = readline(sys.argv[1])
+    dna = buildDNA(line)
+    print dna.printPDB()
+    
+if __name__ == "__main__":
+    main()

DNA/computeParametersOfDNA.py

+from ptools import *
+import sys
+
+nargs = len(sys.argv)
+    if nargs < 2:
+        print "usage: computeParametersOfDNA.py DNA.pdb "
+        raise SystemExit
+
+dna = DNA ("pb.aa.pdb",sys.argv[1])
+
+
+n= len(dna.printParam().split("\n"))
+
+twist=0
+for i in dna.printParam().split("\n"):
+	twist +=float(i.split()[5])
+
+
+roll=0
+for i in dna.printParam().split("\n"):
+	roll +=float(i.split()[8])
+
+
+
+tilt=0
+for i in dna.printParam().split("\n"):
+
+	tilt +=float(i.split()[11])
+
+
+rise=0
+for i in dna.printParam().split("\n"):
+
+	rise +=float(i.split()[14])
+	
+	
+slide=0
+for i in dna.printParam().split("\n"):
+
+	slide +=float(i.split()[17])
+	
+shift=0
+for i in dna.printParam().split("\n"):
+
+	shift +=float(i.split()[20])
+print "Twist(",twist/n,")+Roll(",roll/n,")+Tilt(",tilt/n,")+Rise(",rise/n,")+Slide(",slide/n,")+Shift(",shift/n,")"
+ATOM     67  C1'   A A   3      -0.688   5.757  -0.021 1.60 0.295
+ATOM     68  C2'   A A   3      -0.564   6.922   0.955 1.60-0.067
+ATOM     69  C3'   A A   3      -1.375   8.048   0.314 1.60 0.046
+ATOM     70  C4'   A A   3      -1.116   7.784  -1.168 1.60 0.126
+ATOM     71  O4'   A A   3      -0.896   6.336  -1.300 1.40-0.334
+ATOM     72  C5'   A A   3       0.038   8.551  -1.780 1.60-0.014
+ATOM     73  H1'   A A   3      -1.586   5.183   0.208 1.20 0.064
+ATOM     74  H2'   A A   3       0.466   7.257   1.072 1.20 0.052
+ATOM     75  H2"   A A   3      -1.061   6.606   1.872 1.20 0.052
+ATOM     76  H3'   A A   3      -0.993   9.033   0.583 1.20 0.063
+ATOM     77  H4'   A A   3      -2.026   7.986  -1.733 1.20 0.060
+ATOM     78  H5'   A A   3      -0.282   9.565  -2.020 1.20 0.067
+ATOM     79  H5"   A A   3       0.370   8.047  -2.687 1.20 0.067
+ATOM     80  O3'   A A   3      -2.762   7.934   0.609 1.40-0.281
+ATOM     81  P     A A   3      -3.355   8.095   2.087 1.90 0.456
+ATOM     82  O5'   A A   3      -3.705   6.550   2.306 1.40-0.278
+ATOM     83  O1P   A A   3      -2.298   8.528   3.029 1.40-0.330
+ATOM     84  O2P   A A   3      -4.587   8.914   2.070 1.40-0.330
+ATOM     85  N9    A A   3       0.493   4.850  -0.028 1.50-0.189
+ATOM     86  C8    A A   3       1.832   5.168  -0.058 1.60 0.263
+ATOM     87  N3    A A   3      -0.647   2.703   0.029 1.50-0.411
+ATOM     88  C4    A A   3       0.455   3.480  -0.004 1.60 0.226
+ATOM     89  C5    A A   3       1.753   3.054  -0.021 1.60 0.033
+ATOM     90  C6    A A   3       1.971   1.667  -0.002 1.60 0.361
+ATOM     91  N6    A A   3       3.190   1.109  -0.015 1.50-0.565
+ATOM     92  N7    A A   3       2.619   4.137  -0.055 1.50-0.457
+ATOM     93  C2    A A   3      -0.322   1.430   0.044 1.60 0.249
+ATOM     94  N1    A A   3       0.886   0.870   0.032 1.50-0.451
+ATOM     95  HC2   A A   3      -1.160   0.733   0.071 1.20 0.068
+ATOM     96  H1N   A A   3       3.287   0.094   0.000 1.20 0.295
+ATOM     97  H2N   A A   3       4.021   1.700  -0.041 1.20 0.295
+ATOM     98  HC8   A A   3       2.196   6.195  -0.082 1.20 0.067
+ATOM    695  C1'   T A  23      -0.688  -5.107  -0.021 1.60 0.295
+ATOM    696  C2'   T A  23      -0.539  -6.259  -1.010 1.60-0.067
+ATOM    697  C3'   T A  23      -1.426  -7.352  -0.416 1.60 0.046
+ATOM    698  C4'   T A  23      -1.201  -7.142   1.081 1.60 0.126
+ATOM    699  O4'   T A  23      -0.927  -5.711   1.242 1.40-0.334
+ATOM    700  C5'   T A  23      -0.095  -7.970   1.704 1.60-0.014
+ATOM    701  H1'   T A  23      -1.580  -4.530  -0.262 1.20 0.064
+ATOM    702  H2'   T A  23       0.485  -6.627  -1.073 1.20 0.052
+ATOM    703  H2"   T A  23      -0.973  -5.918  -1.950 1.20 0.052
+ATOM    704  H3'   T A  23      -1.092  -8.349  -0.699 1.20 0.063
+ATOM    705  H4'   T A  23      -2.133  -7.318   1.618 1.20 0.060
+ATOM    706  H5'   T A  23      -0.462  -8.978   1.901 1.20 0.067
+ATOM    707  H5"   T A  23       0.221  -7.507   2.639 1.20 0.067
+ATOM    708  O3'   T A  23      -2.797  -7.154  -0.745 1.40-0.281
+ATOM    709  P     T A  23      -3.382  -7.375  -2.219 1.90 0.456
+ATOM    710  O5'   T A  23      -3.801  -5.854  -2.489 1.40-0.278
+ATOM    711  O1P   T A  23      -2.306  -7.789  -3.147 1.40-0.330
+ATOM    712  O2P   T A  23      -4.577  -8.244  -2.178 1.40-0.330
+ATOM    713  N1    T A  23       0.492  -4.200   0.031 1.50-0.278
+ATOM    714  C6    T A  23       1.763  -4.704   0.111 1.60 0.084
+ATOM    715  O2    T A  23      -0.881  -2.384  -0.077 1.40-0.445
+ATOM    716  N3    T A  23       1.360  -2.038   0.044 1.50-0.404
+ATOM    717  C4    T A  23       2.673  -2.457   0.125 1.60 0.427
+ATOM    718  O4    T A  23       3.592  -1.636   0.163 1.40-0.419
+ATOM    719  C5    T A  23       2.834  -3.893   0.158 1.60-0.059
+ATOM    720  C2    T A  23       0.244  -2.849  -0.005 1.60 0.522
+ATOM    721  C7    T A  23       4.229  -4.433   0.246 1.60-0.147
+ATOM    722  HN3   T A  23       1.200  -1.031   0.019 1.20 0.293
+ATOM    723  HC6   T A  23       1.911  -5.783   0.138 1.20 0.059
+ATOM    724  H1C   T A  23       4.226  -5.491  -0.018 1.20 0.051
+ATOM    725  H2C   T A  23       4.874  -3.888  -0.444 1.20 0.051
+ATOM    726  H3C   T A  23       4.601  -4.314   1.264 1.20 0.051
+ATOM    129  C1'   T T   5      -0.702   5.100  -0.021 1.60 0.295
+ATOM    130  C2'   T T   5      -0.608   6.244   0.982 1.60-0.067
+ATOM    131  C3'   T T   5      -1.395   7.388   0.343 1.60 0.046
+ATOM    132  C4'   T T   5      -1.117   7.149  -1.141 1.60 0.126
+ATOM    133  O4'   T T   5      -0.870   5.710  -1.291 1.40-0.334
+ATOM    134  C5'   T T   5       0.031   7.945  -1.729 1.60-0.014
+ATOM    135  H1'   T T   5      -1.606   4.521   0.167 1.20 0.064
+ATOM    136  H2'   T T   5       0.420   6.574   1.137 1.20 0.052
+ATOM    137  H2"   T T   5      -1.133   5.911   1.875 1.20 0.052
+ATOM    138  H3'   T T   5      -1.009   8.364   0.634 1.20 0.063
+ATOM    139  H4'   T T   5      -2.024   7.343  -1.713 1.20 0.060
+ATOM    140  H5'   T T   5      -0.310   8.951  -1.979 1.20 0.067
+ATOM    141  H5"   T T   5       0.397   7.448  -2.627 1.20 0.067
+ATOM    142  O3'   T T   5      -2.789   7.282   0.613 1.40-0.281
+ATOM    143  P     T T   5      -3.403   7.443   2.081 1.90 0.456
+ATOM    144  O5'   T T   5      -3.761   5.897   2.288 1.40-0.278
+ATOM    145  O1P   T T   5      -2.358   7.871   3.040 1.40-0.330
+ATOM    146  O2P   T T   5      -4.631   8.265   2.049 1.40-0.330
+ATOM    147  N1    T T   5       0.480   4.192  -0.013 1.50-0.278
+ATOM    148  C6    T T   5       1.753   4.695  -0.028 1.60 0.084
+ATOM    149  O2    T T   5      -0.897   2.376   0.025 1.40-0.445
+ATOM    150  N3    T T   5       1.347   2.030   0.018 1.50-0.404
+ATOM    151  C4    T T   5       2.663   2.449   0.004 1.60 0.427
+ATOM    152  O4    T T   5       3.582   1.629   0.012 1.40-0.419
+ATOM    153  C5    T T   5       2.824   3.884  -0.021 1.60-0.059
+ATOM    154  C2    T T   5       0.230   2.841   0.011 1.60 0.522
+ATOM    155  C7    T T   5       4.223   4.425  -0.038 1.60-0.147
+ATOM    156  HN3   T T   5       1.187   1.023   0.036 1.20 0.293
+ATOM    157  HC6   T T   5       1.901   5.775  -0.048 1.20 0.059
+ATOM    158  H1C   T T   5       4.202   5.492   0.187 1.20 0.051
+ATOM    159  H2C   T T   5       4.822   3.908   0.711 1.20 0.051
+ATOM    160  H3C   T T   5       4.660   4.270  -1.024 1.20 0.051
+ATOM    760  C1'   A T  25      -0.702  -5.765  -0.021 1.60 0.295
+ATOM    761  C2'   A T  25      -0.523  -6.941  -0.975 1.60-0.067
+ATOM    762  C3'   A T  25      -1.377  -8.057  -0.375 1.60 0.046
+ATOM    763  C4'   A T  25      -1.195  -7.783   1.119 1.60 0.126
+ATOM    764  O4'   A T  25      -0.984  -6.335   1.248 1.40-0.334
+ATOM    765  C5'   A T  25      -0.073  -8.547   1.792 1.60-0.014
+ATOM    766  H1'   A T  25      -1.585  -5.193  -0.306 1.20 0.064
+ATOM    767  H2'   A T  25       0.512  -7.280  -1.024 1.20 0.052
+ATOM    768  H2"   A T  25      -0.960  -6.633  -1.925 1.20 0.052
+ATOM    769  H3'   A T  25      -0.990  -9.047  -0.617 1.20 0.063
+ATOM    770  H4'   A T  25      -2.133  -7.982   1.637 1.20 0.060
+ATOM    771  H5'   A T  25      -0.408  -9.554   2.039 1.20 0.067
+ATOM    772  H5"   A T  25       0.226  -8.026   2.703 1.20 0.067
+ATOM    773  O3'   A T  25      -2.747  -7.932  -0.745 1.40-0.281
+ATOM    774  P     A T  25      -3.265  -8.082  -2.251 1.90 0.456
+ATOM    775  O5'   A T  25      -3.585  -6.530  -2.475 1.40-0.278
+ATOM    776  O1P   A T  25      -2.168  -8.523  -3.141 1.40-0.330
+ATOM    777  O2P   A T  25      -4.506  -8.885  -2.301 1.40-0.330
+ATOM    778  N9    A T  25       0.478  -4.857   0.046 1.50-0.189
+ATOM    779  C8    A T  25       1.813  -5.176   0.144 1.60 0.263
+ATOM    780  N3    A T  25      -0.658  -2.711  -0.069 1.50-0.411
+ATOM    781  C4    A T  25       0.441  -3.488   0.021 1.60 0.226
+ATOM    782  C5    A T  25       1.736  -3.062   0.103 1.60 0.033
+ATOM    783  C6    A T  25       1.955  -1.675   0.095 1.60 0.361
+ATOM    784  N6    A T  25       3.171  -1.117   0.170 1.50-0.565
+ATOM    785  N7    A T  25       2.599  -4.145   0.181 1.50-0.457
+ATOM    786  C2    A T  25      -0.333  -1.438  -0.068 1.60 0.249
+ATOM    787  N1    A T  25       0.874  -0.878   0.006 1.50-0.451
+ATOM    788  HC2   A T  25      -1.168  -0.741  -0.137 1.20 0.068
+ATOM    789  H1N   A T  25       3.270  -0.102   0.160 1.20 0.295
+ATOM    790  H2N   A T  25       4.000  -1.708   0.238 1.20 0.295
+ATOM    791  HC8   A T  25       2.175  -6.203   0.186 1.20 0.067
+ATOM     99  C1'   C C   4      -0.410   5.271  -0.003 1.60 0.295
+ATOM    100  C2'   C C   4      -0.312   6.430   0.983 1.60-0.067
+ATOM    101  C3'   C C   4      -1.107   7.569   0.346 1.60 0.046
+ATOM    102  C4'   C C   4      -0.832   7.321  -1.137 1.60 0.126
+ATOM    103  O4'   C C   4      -0.586   5.881  -1.272 1.40-0.334
+ATOM    104  C5'   C C   4       0.314   8.113  -1.735 1.60-0.014
+ATOM    105  H1'   C C   4      -1.313   4.694   0.194 1.20 0.064
+ATOM    106  H2'   C C   4       0.716   6.762   1.126 1.20 0.052
+ATOM    107  H2"   C C   4      -0.828   6.104   1.886 1.20 0.052
+ATOM    108  H3'   C C   4      -0.723   8.548   0.630 1.20 0.063
+ATOM    109  H4'   C C   4      -1.741   7.508  -1.709 1.20 0.060
+ATOM    110  H5'   C C   4      -0.023   9.121  -1.976 1.20 0.067
+ATOM    111  H5"   C C   4       0.668   7.618  -2.639 1.20 0.067
+ATOM    112  O3'   C C   4      -2.499   7.459   0.623 1.40-0.281
+ATOM    113  P     C C   4      -3.105   7.622   2.094 1.90 0.456
+ATOM    114  O5'   C C   4      -3.468   6.079   2.316 1.40-0.278
+ATOM    115  O1P   C C   4      -2.056   8.051   3.045 1.40-0.330
+ATOM    116  O2P   C C   4      -4.332   8.446   2.065 1.40-0.330
+ATOM    117  N1    C C   4       0.771   4.363   0.005 1.50-0.240
+ATOM    118  C6    C C   4       2.037   4.859  -0.011 1.60 0.252
+ATOM    119  O2    C C   4      -0.635   2.587   0.043 1.40-0.447
+ATOM    120  N3    C C   4       1.598   2.147   0.037 1.50-0.518
+ATOM    121  C4    C C   4       2.848   2.624   0.022 1.60 0.407
+ATOM    122  N4    C C   4       3.851   1.759   0.030 1.50-0.560
+ATOM    123  C5    C C   4       3.113   4.031  -0.003 1.60-0.279
+ATOM    124  C2    C C   4       0.534   2.991   0.029 1.60 0.469
+ATOM    125  H1N   C C   4       3.660   0.757   0.048 1.20 0.295
+ATOM    126  H2N   C C   4       4.814   2.095   0.020 1.20 0.295
+ATOM    127  HC5   C C   4       4.132   4.421  -0.015 1.20 0.055
+ATOM    128  HC6   C C   4       2.191   5.938  -0.030 1.20 0.060
+ATOM    727  C1'   G C  24      -0.410  -5.594  -0.003 1.60 0.295
+ATOM    728  C2'   G C  24      -0.238  -6.769  -0.959 1.60-0.067
+ATOM    729  C3'   G C  24      -1.102  -7.884  -0.373 1.60 0.046
+ATOM    730  C4'   G C  24      -0.918  -7.623   1.123 1.60 0.126
+ATOM    731  O4'   G C  24      -0.684  -6.180   1.260 1.40-0.334
+ATOM    732  C5'   G C  24       0.192  -8.407   1.793 1.60-0.014
+ATOM    733  H1'   G C  24      -1.294  -5.020  -0.279 1.20 0.064
+ATOM    734  H2'   G C  24       0.795  -7.116  -1.005 1.20 0.052
+ATOM    735  H2"   G C  24      -0.668  -6.453  -1.909 1.20 0.052
+ATOM    736  H3'   G C  24      -0.723  -8.874  -0.624 1.20 0.063
+ATOM    737  H4'   G C  24      -1.859  -7.811   1.639 1.20 0.060
+ATOM    738  H5'   G C  24      -0.154  -9.416   2.019 1.20 0.067
+ATOM    739  H5"   G C  24       0.487  -7.906   2.714 1.20 0.067
+ATOM    740  O3'   G C  24      -2.469  -7.746  -0.743 1.40-0.281
+ATOM    741  P     G C  24      -2.985  -7.918  -2.248 1.90 0.456
+ATOM    742  O5'   G C  24      -3.351  -6.383  -2.510 1.40-0.278
+ATOM    743  O1P   G C  24      -1.874  -8.343  -3.131 1.40-0.330
+ATOM    744  O2P   G C  24      -4.203  -8.756  -2.287 1.40-0.330
+ATOM    745  N9    G C  24       0.769  -4.686   0.064 1.50-0.185
+ATOM    746  C8    G C  24       2.106  -5.001   0.161 1.60 0.227
+ATOM    747  N2    G C  24      -1.157  -0.376  -0.141 1.50-0.561
+ATOM    748  N3    G C  24      -0.412  -2.567  -0.053 1.50-0.533
+ATOM    749  C4    G C  24       0.723  -3.308   0.038 1.60 0.226
+ATOM    750  C5    G C  24       2.025  -2.875   0.120 1.60-0.035
+ATOM    751  C6    G C  24       2.295  -1.481   0.115 1.60 0.432
+ATOM    752  O6    G C  24       3.379  -0.907   0.181 1.40-0.430
+ATOM    753  N7    G C  24       2.899  -3.958   0.199 1.50-0.450
+ATOM    754  C2    G C  24      -0.160  -1.260  -0.058 1.60 0.499
+ATOM    755  N1    G C  24       1.113  -0.732   0.021 1.50-0.360
+ATOM    756  HN1   G C  24       1.202   0.284   0.011 1.20 0.294
+ATOM    757  H1N   G C  24      -0.953   0.624  -0.145 1.20 0.297
+ATOM    758  H2N   G C  24      -2.122  -0.699  -0.203 1.20 0.297
+ATOM    759  HC8   G C  24       2.473  -6.027   0.204 1.20 0.067
+ATOM     34  C1'   G G   2      -0.379   5.591   0.001 1.60 0.295
+ATOM     35  C2'   G G   2      -0.240   6.782   0.944 1.60-0.067
+ATOM     36  C3'   G G   2      -1.062   7.893   0.293 1.60 0.046
+ATOM     37  C4'   G G   2      -0.834   7.597  -1.190 1.60 0.126
+ATOM     38  O4'   G G   2      -0.610   6.146  -1.285 1.40-0.334
+ATOM     39  C5'   G G   2       0.304   8.355  -1.844 1.60-0.014
+ATOM     40  H1'   G G   2      -1.274   5.022   0.257 1.20 0.064
+ATOM     41  H2'   G G   2       0.793   7.118   1.034 1.20 0.052
+ATOM     42  H2"   G G   2      -0.721   6.488   1.877 1.20 0.052
+ATOM     43  H3'   G G   2      -0.673   8.884   0.531 1.20 0.063
+ATOM     44  H4'   G G   2      -1.757   7.781  -1.740 1.20 0.060
+ATOM     45  H5'   G G   2      -0.036   9.356  -2.112 1.20 0.067
+ATOM     46  H5"   G G   2       0.618   7.823  -2.742 1.20 0.067
+ATOM     47  O3'   G G   2      -2.443   7.788   0.617 1.40-0.281
+ATOM     48  P     G G   2      -3.016   7.911   2.107 1.90 0.456
+ATOM     49  O5'   G G   2      -3.331   6.356   2.307 1.40-0.278
+ATOM     50  O1P   G G   2      -1.953   8.352   3.038 1.40-0.330
+ATOM     51  O2P   G G   2      -4.264   8.704   2.122 1.40-0.330
+ATOM     52  N9    G G   2       0.802   4.683  -0.006 1.50-0.185
+ATOM     53  C8    G G   2       2.142   4.999  -0.036 1.60 0.227
+ATOM     54  N2    G G   2      -1.133   0.373   0.102 1.50-0.561
+ATOM     55  N3    G G   2      -0.384   2.565   0.051 1.50-0.533
+ATOM     56  C4    G G   2       0.754   3.306   0.018 1.60 0.226
+ATOM     57  C5    G G   2       2.058   2.873   0.001 1.60-0.035
+ATOM     58  C6    G G   2       2.328   1.479   0.019 1.60 0.432
+ATOM     59  O6    G G   2       3.413   0.905   0.007 1.40-0.430
+ATOM     60  N7    G G   2       2.935   3.955  -0.035 1.50-0.450
+ATOM     61  C2    G G   2      -0.133   1.257   0.068 1.60 0.499
+ATOM     62  N1    G G   2       1.143   0.730   0.053 1.50-0.360
+ATOM     63  HN1   G G   2       1.231  -0.286   0.068 1.20 0.294
+ATOM     64  H1N   G G   2      -0.929  -0.627   0.115 1.20 0.297
+ATOM     65  H2N   G G   2      -2.100   0.696   0.114 1.20 0.297
+ATOM     66  HC8   G G   2       2.510   6.024  -0.061 1.20 0.067
+ATOM    665  C1'   C G  22      -0.379  -5.274   0.001 1.60 0.295
+ATOM    666  C2'   C G  22      -0.233  -6.432  -0.980 1.60-0.067
+ATOM    667  C3'   C G  22      -1.206  -7.476  -0.436 1.60 0.046
+ATOM    668  C4'   C G  22      -0.987  -7.296   1.066 1.60 0.126
+ATOM    669  O4'   C G  22      -0.618  -5.886   1.260 1.40-0.334
+ATOM    670  C5'   C G  22       0.050  -8.206   1.692 1.60-0.014
+ATOM    671  H1'   C G  22      -1.272  -4.696  -0.240 1.20 0.064
+ATOM    672  H2'   C G  22       0.773  -6.851  -0.987 1.20 0.052
+ATOM    673  H2"   C G  22      -0.596  -6.071  -1.942 1.20 0.052
+ATOM    674  H3'   C G  22      -0.932  -8.489  -0.731 1.20 0.063
+ATOM    675  H4'   C G  22      -1.938  -7.416   1.586 1.20 0.060
+ATOM    676  H5'   C G  22      -0.360  -9.208   1.812 1.20 0.067
+ATOM    677  H5"   C G  22       0.345  -7.812   2.664 1.20 0.067
+ATOM    678  O3'   C G  22      -2.554  -7.182  -0.788 1.40-0.281
+ATOM    679  P     C G  22      -3.135  -7.518  -2.241 1.90 0.456
+ATOM    680  O5'   C G  22      -3.650  -6.066  -2.669 1.40-0.278
+ATOM    681  O1P   C G  22      -2.036  -7.947  -3.137 1.40-0.330
+ATOM    682  O2P   C G  22      -4.280  -8.447  -2.126 1.40-0.330
+ATOM    683  N1    C G  22       0.801  -4.366   0.052 1.50-0.240
+ATOM    684  C6    C G  22       2.064  -4.862   0.131 1.60 0.252
+ATOM    685  O2    C G  22      -0.602  -2.589  -0.056 1.40-0.447
+ATOM    686  N3    C G  22       1.627  -2.149   0.062 1.50-0.518
+ATOM    687  C4    C G  22       2.875  -2.626   0.140 1.60 0.407
+ATOM    688  N4    C G  22       3.877  -1.762   0.182 1.50-0.560
+ATOM    689  C5    C G  22       3.139  -4.035   0.178 1.60-0.279
+ATOM    690  C2    C G  22       0.565  -2.994   0.016 1.60 0.469
+ATOM    691  H1N   C G  22       3.688  -0.761   0.155 1.20 0.295
+ATOM    692  H2N   C G  22       4.839  -2.098   0.242 1.20 0.295
+ATOM    693  HC5   C G  22       4.155  -4.424   0.242 1.20 0.055
+ATOM    694  HC6   C G  22       2.216  -5.941   0.158 1.20 0.060
+ATOM      1  GP1   A A   3      -2.982   8.060   2.598   30  -0.700 0 0
+ATOM      2  GS1   A A   3      -0.931   7.674  -0.352   31   0.000 0 0
+ATOM      3  GS2   A A   3      -0.337   5.927  -0.008   32   0.000 0 0
+ATOM      4  GA1   A A   3       2.028   4.306  -0.008   33   0.000 0 0
+ATOM      5  GA2   A A   3       2.330   1.263  -0.017   34   0.000 0 0
+ATOM      6  GA3   A A   3       0.210   2.125   0.064   35   0.000 0 0
+ATOM      7  GP1   T A  23      -3.094  -8.174  -2.570   30  -0.700 0 0
+ATOM      8  GS1   T A  23      -0.992  -7.823   0.314   31   0.000 0 0
+ATOM      9  GS2   T A  23      -0.337  -6.105  -0.008   32   0.000 0 0
+ATOM     10  GT1   T A  23       0.784  -3.946  -0.052   41   0.000 0 0
+ATOM     11  GT2   T A  23       3.320  -3.306   0.039   42   0.000 0 0
+ATOM      1  GP1   T T   5      -3.017   8.225   2.614   30  -0.700 0 0
+ATOM      2  GS1   T T   5      -0.951   7.855  -0.321   31   0.000 0 0
+ATOM      3  GS2   T T   5      -0.349   6.093   0.003   32   0.000 0 0
+ATOM      4  GT1   T T   5       0.772   3.934   0.046   41   0.000 0 0
+ATOM      5  GT2   T T   5       3.308   3.294  -0.039   42   0.000 0 0
+ATOM      6  GP1   A T  25      -2.998  -8.049  -2.606   30  -0.700 0 0
+ATOM      7  GS1   A T  25      -0.950  -7.679   0.356   31   0.000 0 0
+ATOM      8  GS2   A T  25      -0.349  -5.940   0.003   32   0.000 0 0
+ATOM      9  GA1   A T  25       2.017  -4.319   0.003   33   0.000 0 0
+ATOM     10  GA2   A T  25       2.319  -1.275   0.014   34   0.000 0 0
+ATOM     11  GA3   A T  25       0.199  -2.138  -0.073   35   0.000 0 0
+ATOM      1  GP1   C C   4      -2.935   8.225   2.616   30  -0.700 0 0
+ATOM      2  GS1   C C   4      -0.886   7.849  -0.328   31   0.000 0 0
+ATOM      3  GS2   C C   4      -0.280   6.090   0.005   32   0.000 0 0
+ATOM      4  GC1   C C   4       0.837   3.928   0.052   39   0.000 0 0
+ATOM      5  GC2   C C   4       3.012   2.869  -0.023   40   0.000 0 0
+ATOM      6  GP1   G C  24      -2.932  -8.070  -2.609   30  -0.700 0 0
+ATOM      7  GS1   G C  24      -0.900  -7.694   0.345   31   0.000 0 0
+ATOM      8  GS2   G C  24      -0.280  -5.954   0.005   32   0.000 0 0
+ATOM      9  GG1   G C  24       2.089  -4.321   0.005   36   0.000 0 0
+ATOM     10  GG2   G C  24       2.366  -1.271   0.016   37   0.000 0 0
+ATOM     11  GG3   G C  24      -0.093  -1.649  -0.083   38   0.000 0 0
+ATOM      1  GP1   G G   2      -2.861   8.043   2.635   30  -0.700 0 0
+ATOM      2  GS1   G G   2      -0.867   7.680  -0.376   31   0.000 0 0
+ATOM      3  GS2   G G   2      -0.258   5.938  -0.003   32   0.000 0 0
+ATOM      4  GG1   G G   2       2.112   4.309  -0.003   36   0.000 0 0
+ATOM      5  GG2   G G   2       2.388   1.260  -0.016   37   0.000 0 0
+ATOM      6  GG3   G G   2      -0.070   1.638   0.091   38   0.000 0 0
+ATOM      7  GP1   C G  22      -3.063  -8.164  -2.565   30  -0.700 0 0
+ATOM      8  GS1   C G  22      -1.015  -7.782   0.276   31   0.000 0 0
+ATOM      9  GS2   C G  22      -0.258  -6.101  -0.003   32   0.000 0 0
+ATOM     10  GC1   C G  22       0.858  -3.939  -0.052   39   0.000 0 0
+ATOM     11  GC2   C G  22       3.034  -2.881   0.016   40   0.000 0 0

DNA/rotateDNAaroundAxis.py

+from ptools import *
+import math
+import sys
+import os
+
+PB_CG = "pb.cg.pdb"
+PB_AA = "pb.aa.pdb"
+
+def ScalProd(a,b):
+    return a.x * b.x + a.y * b.y + a.z * b.z ;
+
+def readline(lineName):
+    line=[]
+    l = Rigidbody(lineName)
+    for i in xrange(l.Size()):
+        line.append(l.CopyAtom(i))
+    return line
+   
+def nextModel(line,D,pos):
+    start = line [pos]
+    nextpos = None
+    for i,dot in enumerate(line[pos+1:]):
+        if Dist (dot,start)>D:
+            nextpos=pos+1+i-1
+            A= line[pos+1+i-1]
+            B= line[pos+1+i]
+            break
+    if not nextpos:
+        return None,pos
+    AB = B.GetCoords () - A.GetCoords()
+    Astart = start.GetCoords () - A.GetCoords()
+    AB.Normalize()
+    dproj = ScalProd (AB,Astart)
+    proj = Atom(Atomproperty (),A.GetCoords()+(AB*dproj))
+    startProj = Dist(start,proj)
+    distX = math.sqrt((D**2)-(startProj**2))#V
+    at = Atomproperty ()
+    at.SetAtomId (pos)
+    X= Atom(at,proj.GetCoords()+(AB*distX))
+
+    
+    model =Rigidbody()
+    model.AddAtom ( start )
+    model.AddAtom ( X )
+    return model,nextpos
+        
+def getAxislist(line):
+    segment = DNA(PB_AA ,"AA",BDNA())
+    center0 =Atom(Atomproperty (), segment[0].getRigidBody ().FindCenter ())
+    center1 =Atom(Atomproperty (), segment[1].getRigidBody ().FindCenter ())
+    D = Dist(center0,center1)
+    pos = 0
+    
+    #extract the axis of rotation from the line
+    axislist = []
+    model,newpos = nextModel(line,D,pos)
+    axislist.append(model)
+    while True :
+        pos = newpos
+        model,newpos = nextModel(line,D,pos)
+        if newpos == pos:
+            break
+        print model.PrintPDB()
+        axislist.append(model)
+    return axislist
+     
+def rotateDNA(dna,axislist,angle):
+    new = AttractRigidbody()
+    for i,model in enumerate(axislist):
+        bp = dna[i].getRigidBody()
+        
+        #rotation
+        ABrotate( model.CopyAtom(0).GetCoords(), model.CopyAtom(1).GetCoords(), bp, math.radians(angle) )
+        #print bp.PrintPDB()
+        new = new + bp
+        #print new.size()
+    bp = dna[dna.size()-1].getRigidBody()
+    ABrotate( model.CopyAtom(0).GetCoords(), model.CopyAtom(1).GetCoords(), bp, math.radians(angle) )
+    new = new + bp
+    return new
+    #python screwOptimisation.py ../1QRV_93.axis ../1QRV_93.candidate ../../cgPyAttract/1QRV_bound-prot.red
+##python screwOptimisation.py ../3CRO_93.dots.3clust.axis ../3CRO_93.dots.3clust.candidate.pdb ../../cgPyAttract/3CRO_bound-prot.red >toto.pdb                           python screwOptimisation.py ../3CRO_93.dots.3clust.axis 3CRO_test.pdb ../../cgPyAttract/3CRO_bound-prot.red >toto.pdb   
+        
+
+def main():
+    nargs = len(sys.argv)
+    if nargs < 3:
+        print "usage: makeCandidate.py axis DNA.pdb angle(degree)"
+        raise SystemExit
+
+    line = readline(sys.argv[1])
+    dna = DNA(PB_AA, Rigidbody(sys.argv[2]))
+    angle = float(sys.argv[3])
+    
+    axislist = getAxislist(line)
+    
+    d=rotateDNA(dna,axislist,angle)
+    print d.PrintPDB()
+    
+if __name__ == "__main__":
+    main()

Heligeom/buildProteinAlongAnAxis.py

+from ptools import *
+import math
+import sys
+import os
+import string
+
+def readline(lineName):
+    line=[]
+    l = Rigidbody(lineName)
+    for i in xrange(l.Size()):
+        line.append(l.CopyAtom(i))
+    return line
+    
+def changeChain(rig,letter):
+    rsize = rig.Size()
+    for i in xrange(0,rsize):
+        at=rig.GetAtomProperty (i)
+        at.SetChainId(letter)
+        rig.SetAtomProperty (i,at)
+    return rig
+    
+
+def cleanPDB(rig):
+    rSize = rig.Size()
+    l=[]
+    base =[]
+    nb=-1
+    oldChain = None
+    for i in xrange(0,rSize):
+        a=rig.CopyAtom(i)
+        chain = a.GetChainId ()
+        resid = a.GetResidId ()
+        if nb != resid:
+            if nb != -1:
+                l.append([oldChain,base])
+            nb = resid
+            base=[]
+        if oldChain != chain:
+            oldChain = chain
+        
+        base.append(a)
+    l.append([chain,base])
+    
+    l = sorted(l, key=lambda base: base[0])
+    
+    for i,b in enumerate(l):
+        if b[0] == "B":
+            break
+    l[i],l[i+1]=l[i+1],l[i]
+    
+    newRig = Rigidbody()
+    nbAtom=0
+    for nbRes,base in enumerate(l):
+        #base = base[1]
+        
+        for a in base[1]:
+            a.SetResidId (nbRes)
+            a.SetAtomId (nbAtom)
+            nbAtom += 1
+            newRig.AddAtom(a)
+            #print a.ToPdbString(),
+    return newRig
+
+
+def ScalProd(a,b):
+    return a.x * b.x + a.y * b.y + a.z * b.z ;
+
+def VectProd(u,v):
+    UvectV = Coord3D()
+    UvectV.x = u.y * v.z - u.z * v.y 
+    UvectV.y = u.z * v.x - u.x * v.z 
+    UvectV.z = u.x * v.y - u.y * v.x
+    return UvectV 
+
+def nextModel(line,D,pos):
+    start = line [pos]
+    nextpos = None
+    for i,dot in enumerate(line[pos+1:]):
+        if Dist (dot,start)>D:
+            nextpos=pos+1+i-1
+            A= line[pos+1+i-1]
+            B=line[pos+1+i]
+            break
+    if not nextpos:
+        return None,pos
+    AB = B.GetCoords () - A.GetCoords()
+    Astart = start.GetCoords () - A.GetCoords()
+    AB.Normalize()
+    dproj = ScalProd (AB,Astart)
+    proj = Atom(Atomproperty (),A.GetCoords()+(AB*dproj))
+    startProj = Dist(start,proj)
+    distX = math.sqrt((D**2)-(startProj**2))#V
+    at = Atomproperty ()
+    at.SetAtomId (pos)
+    X= Atom(at,proj.GetCoords()+(AB*distX))
+
+    
+    model =Rigidbody()
+    model.AddAtom ( start )
+    model.AddAtom ( X )
+    return model,nextpos
+        
+def buildProt(line,mono1,hp, angle):
+    #place the first monomer on the correct the axis
+    p = hp.point
+    if hp.normtranslation < 0 :
+        v = -1*hp.unitVector.Normalize()
+    else:
+        v = hp.unitVector.Normalize()
+        #line.reverse()
+    D = abs(hp.normtranslation)
+    hpangle = hp.angle
+    #D = (hp.normtranslation)
+    pos = 0
+    model,newpos = nextModel(line,D,pos)
+    
+    
+
+    #build a 3 point ref for the monomer
+    
+    monoCenter = mono1.FindCenter ()
+
+    #find vector
+    Pcenter =  monoCenter - p
+    
+    
+    #scal prod to find projection of mono center on the axis
+    scal = ScalProd(v,Pcenter)
+        
+    #for i in xrange (0,100):
+        #print Atom(Atomproperty(),( p + v*i )).ToPdbString ()
+    refMobil = Rigidbody()
+    O= p+(v*scal)
+    at=Atomproperty()
+    at.SetType ("O")
+    refMobil.AddAtom ( at,O)
+    at.SetType ("Z")
+    refMobil.AddAtom ( at,O+(v))            
+    at.SetType ("Y")
+    refMobil.AddAtom ( at,O+(monoCenter-O).Normalize())
+    v_third = VectProd(monoCenter - O, v)
+    at.SetType ("X")
+    refMobil.AddAtom ( Atom(at,O+v_third.Normalize()))
+    distCenterToAxis = Dist (Atom(at,O), Atom(at,monoCenter))
+    
+    A,B=model.CopyAtom(0).GetCoords (),model.CopyAtom(1).GetCoords ()
+    v = (A-B).Normalize ()
+    
+    refMobil2 = refMobil.SelectAllAtoms().CreateRigid()
+    refMobil2.ABrotate ( A, A + v, hpangle )
+    refMobil2.Translate( v * D)
+    segmentMobil = refMobil + refMobil2
+    
+    #build a 3 point ref for the axis
+    #get unit vector on the first local axis of line
+
+
+    #create a non parrallel vector
+    v2 = Coord3D(v)
+    v2.x = v2.x+1
+    v2.y = v2.y-1
+    #get an orthogonal vector with cross product
+    vortho = VectProd(v,v2)
+    vortho =vortho.Normalize () 
+    
+    #construct the second rigidbody
+    refModel = Rigidbody()
+    refModel.AddAtom ( Atom(Atomproperty(),A))
+    refModel.AddAtom ( Atom(Atomproperty(),A+v))
+    refModel.AddAtom ( Atom(Atomproperty(),A+vortho))
+    
+    v_third = VectProd((A+vortho*distCenterToAxis)-A,v)
+    refModel.AddAtom ( Atom(Atomproperty(),A+v_third.Normalize()))
+    #rotate vortho
+    #temp = Rigidbody()
+    #temp.AddAtom ( Atom(Atomproperty(),A))
+    #temp.AddAtom ( Atom(Atomproperty(),A+vortho))
+    #temp.ABrotate ( A, A + v,angle )
+    #vorthoA,vorthoB=temp.CopyAtom(0).GetCoords (),temp.CopyAtom(1).GetCoords ()
+    #vortho =(vorthoA-vorthoB).Normalize ()
+    #get the mobil
+    mobil =  Rigidbody()
+    mobil.AddAtom ( Atom(Atomproperty(),A))
+    mobil.AddAtom ( Atom(Atomproperty(),A+v*D))
+    
+    #move the monomer
+    refModel.ABrotate ( A, A + v,angle  )
+    m = superpose(refModel,refMobil).matrix
+    #mono1.ApplyMatrix (m)
+    #print refModel.PrintPDB()
+    #mono1.ABrotate ( A, A + v,angle  )
+    #refModel.ABrotate ( A, A + v,angle  )
+    mono1= changeChain(mono1,"A")
+    #print mono1.PrintPDB()
+    #get the matrix
+    refModel2 = refModel.SelectAllAtoms().CreateRigid()
+    refModel2.ABrotate ( A, A + v, hpangle )
+    refModel2.Translate( v * D)
+    #for i in xrange (0,10):
+        #print Atom(Atomproperty(),( A + v*i*3 )).ToPdbString ()
+        
+        
+    segment = refModel + refModel2
+    segment2 = segment.SelectAllAtoms().CreateRigid()
+    segment2.ABrotate ( A, A + v, hpangle )
+    segment2.Translate( v * D)
+    m_rot = superpose(segment2,segment).matrix
+    
+    #print segment.PrintPDB()
+    #print segment2.PrintPDB()
+    size = 2
+    nbChain=1
+    newpos =0
+    while True :
+        pos = newpos
+        model,newpos = nextModel(line,D,pos)
+        if newpos == pos:
+            break
+                
+        #extract next segment and
+        A,B=model.CopyAtom(0).GetCoords (),model.CopyAtom(1).GetCoords ()
+        vector = (A-B).Normalize ()
+        #vector = (B-A).Normalize ()
+        #newRef = mobil.SelectAllAtoms().CreateRigid()
+        
+        
+        v_temp = VectProd(vector,A+vortho*distCenterToAxis)
+        v_temp = VectProd(vector,v_temp)
+        v_temp = v_temp.Normalize ()
+        newRef =Rigidbody()
+        at=Atomproperty()
+        at.SetChainId (string.ascii_uppercase[nbChain%26])
+        at.SetType ("O")
+        newRef.AddAtom (at,A)
+        at.SetType ("Z")
+        newRef.AddAtom (at,A+vector)
+        at.SetType ("Y")
+        newRef.AddAtom ( Atom(at,A+v_temp))
+        at.SetType ("X")    
+        v_third = VectProd((A+v_temp*distCenterToAxis)-A,vector)
+        newRef.AddAtom ( Atom(at,A+v_third.Normalize()))
+        ##vector = (A-B).Normalize ()
+        #newRef.ABrotate ( A, A + vector, angle )
+        #segment2.Translate( vector * hp.normtranslation )
+         
+        
+        for i in xrange(0,size-1):
+            newRef.ABrotate ( A, A + vector, hpangle )
+        newRef.ABrotate ( A, A + vector, angle )            
+        size+=1
+        
+        #print newRef.PrintPDB()
+        newRef2 = newRef.SelectAllAtoms().CreateRigid()
+        newRef2.ABrotate ( A, A + vector, hpangle )
+        newRef2.Translate( vector * D)
+    
+        #newRef += newRef2
+        #refMobil
+        m = superpose(newRef,refMobil).matrix
+        #print superpose(newRef,refMobil).rmsd
+        #m = superpose(newRef,segment).matrix
+        newMono = mono1.SelectAllAtoms().CreateRigid()
+        #newMono =  BasePair (changeChain(newMono,string.ascii_uppercase[nbChain%26]))
+        newMono = (changeChain(newMono,string.ascii_uppercase[nbChain%26]))
+        
+        #newMono.apply(m)
+        newMono.ApplyMatrix(m)
+        nbChain+=1
+        print newMono.PrintPDB()
+        mob = refMobil.SelectAllAtoms().CreateRigid()
+        mob.ApplyMatrix(m)
+        mob= (changeChain(mob,string.ascii_uppercase[nbChain%26]))
+        #print mob.PrintPDB()
+        
+        #print refModel2.getRigidBody ().PrintPDB()
+    
+def main():
+    nargs = len(sys.argv)
+    if nargs < 4:
+        print "usage: buildProteinAlongAnAxis.py  axis.pdb monomer1.pdb monomer2.pdb [angle degree]"
+        raise SystemExit
+
+    line = readline(sys.argv[1])
+    mono1 = Rigidbody( sys.argv[2])
+    mono2 = Rigidbody( sys.argv[3])
+    
+    angle = math.radians(0)
+    if nargs >4 :
+        angle = math.radians(float(sys.argv[4]))
+    
+    hp = MatTrans2screw(superpose (mono2,mono1).matrix)
+    
+    prot = buildProt(line,mono1,hp,angle)
+    
+    #for i in xrange(0,360,6):
+        #hp = MatTrans2screw(superpose (mono2,mono1).matrix)
+        #mono = mono1.SelectAllAtoms().CreateRigid()
+        #print "MODEL       "+str(i)
+        #prot = buildProt(line,mono,hp,math.radians(float(i)))
+        #print "ENDMDL"  
+    
+if __name__ == "__main__":
+    main()

Heligeom/extractHelicoidalModel.py

+#!/usr/bin/env python