SSAGES  0.8.5
Software Suite for Advanced General Ensemble Simulations
ReadBackbone.h
1 
21 #pragma once
22 
23 #include <fstream>
24 
25 // Not a full PDB parser!!! Only grabs backbone atoms given residue numbers and
26 // a reference PDB. This will not parse any other PDB record type other than
27 // "ATOM", and does not read any info from a PDB file past the column containing
28 // residue number. Could be extended to parse a PDB file completely, which would
29 // be better placed in ReadFile.h, but sticking with minimum requirements for
30 // secondary structure CVs here.
31 
32 namespace SSAGES
33 {
35 
43  {
44  public:
46 
56  //static std::vector<int> GetPdbBackbone(std::string refpdb, std::vector<int> resids)
57  static std::vector< std::vector< std::string> > GetPdbBackbone(std::string refpdb, std::vector<int> resids)
58  {
59  //std::vector<int> atomids;
60  std::vector< std::vector< std::string > > atomids(2, std::vector<std::string> (0));
61  //std::vector<int> atomnums;
62  std::vector<std::string> atomnums;
63  std::vector<std::string> atomtypes;
64  std::vector<std::string> resnames;
65  std::vector<std::string> chains;
66  std::vector<int> resnums;
67  //std::vector<int> tempatoms(5);
68  std::vector<std::string> tempatoms(5);
69  std::vector<std::string> backboneAtoms = {"N", "CA", "CB", "C", "O", "OT1"};
70  std::vector<std::string> glycineAtoms = {"N", "CA", "HA1", "C", "O", "OT1"};
71  std::string atomnum, resnum;
72  std::string record, atomtype, resname, chain;
73  std::ifstream pdbfile;
74  pdbfile.open(refpdb, std::ios::in);
75  std::string line;
76 
77  while(std::getline(pdbfile, line)){
78  if(line.length() < 26) line.append(26 - line.length(), ' ');
79  record = line.substr(0, 6);
80  atomnum = line.substr(6, 5);
81  atomtype = line.substr(12, 4);
82  resname = line.substr(17, 3);
83  chain = line.substr(21, 1);
84  resnum = line.substr(22, 4);
85  record.erase( std::remove( record.begin(), record.end(), ' '), record.end());
86  if((record == "ATOM") && ((std::stoi(resnum) - resids[0]) <= (resids.back() - resids[0]))){
87  atomtype.erase( std::remove( atomtype.begin(), atomtype.end(), ' '), atomtype.end());
88  resname.erase( std::remove( resname.begin(), resname.end(), ' '), resname.end());
89  if(resname == "GLY" && std::find(glycineAtoms.begin(), glycineAtoms.end(), atomtype) != glycineAtoms.end()){
90  //atomnums.push_back(std::stoi(atomnum));
91  atomnums.push_back(atomnum);
92  atomtypes.push_back(atomtype);
93  resnames.push_back(resname);
94  chains.push_back(chain);
95  resnums.push_back(std::stoi(resnum));
96  } else if(std::find(backboneAtoms.begin(), backboneAtoms.end(), atomtype) != backboneAtoms.end()){
97  //atomnums.push_back(std::stoi(atomnum));
98  atomnums.push_back(atomnum);
99  atomtypes.push_back(atomtype);
100  resnames.push_back(resname);
101  chains.push_back(chain);
102  resnums.push_back(std::stoi(resnum));
103  }
104  }
105  }
106 
107  // atoms in PDB not necessarily in N CA CB C O order, fix that:
108  for(size_t i = 0; i < resids.size(); i++){
109  for(size_t j = 5 * i; j < (5 * i + 5); j++){
110  if(atomtypes[j] == "N"){
111  tempatoms[0] = atomnums[j];
112  std::cout << "N - Atom " << tempatoms[0] << " - Chain " << chains[j] << std::endl;
113  } else if(atomtypes[j] == "CA"){
114  tempatoms[1] = atomnums[j];
115  std::cout << "CA - Atom " << tempatoms[1] << " - Chain " << chains[j] << std::endl;
116  } else if(atomtypes[j] == "C"){
117  tempatoms[3] = atomnums[j];
118  std::cout << "C - Atom " << tempatoms[3] << " - Chain " << chains[j] << std::endl;
119  } else if(atomtypes[j] == "O"){
120  tempatoms[4] = atomnums[j];
121  std::cout << "O - Atom " << tempatoms[4] << " - Chain " << chains[j] << std::endl;
122  } else if(atomtypes[j] == "OT1"){
123  tempatoms[4] = atomnums[j];
124  std::cout << "OT1- Atom " << tempatoms[4] << " - Chain " << chains[j] << std::endl;
125  } else{
126  tempatoms[2] = atomnums[j];
127  std::cout << atomtypes[j] << " - Atom " << tempatoms[2] << " - Chain " << chains[j] << std::endl;
128  }
129  atomids[1].push_back(chains[j]);
130  }
131  //atomids.insert(atomids.end(), tempatoms.begin(), tempatoms.end());
132  atomids[0].insert(atomids[0].end(), tempatoms.begin(), tempatoms.end());
133  }
134 
135  return atomids;
136  }
137  };
138 }
static std::vector< std::vector< std::string > > GetPdbBackbone(std::string refpdb, std::vector< int > resids)
Read protein backbone atoms from reference PDB file.
Definition: ReadBackbone.h:57
Utility class to read protein backbone atoms from a reference file.
Definition: ReadBackbone.h:42