diff --git a/Biopool/APPS/Cif2Secondary.cc b/Biopool/APPS/Cif2Secondary.cc new file mode 100644 index 0000000..5db0500 --- /dev/null +++ b/Biopool/APPS/Cif2Secondary.cc @@ -0,0 +1,115 @@ +/* This file is part of Victor. + + Victor is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + Victor 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 General Public License for more details. + + You should have received a copy of the GNU General Public License + along with Victor. If not, see . + */ + +#include +#include +#include +#include +#include + +using namespace Victor; +using namespace Victor::Biopool; + +void sShowHelp() { + cout << "Cif 2 Secondary Structure converter\n" + << "\t H = helix, \t E = extended (strand, sheet), \t . = other.\n" + << " Options: \n" + << "\t-i \t\t Input file for CIF structure\n" + << "\n"; +} + +int main(int argc, char** argv) { + + if (getArg("h", argc, argv)) { + sShowHelp(); + return 1; + }; + vector allCh; + string chainID = "!"; + string inputFile; + getArg("i", inputFile, argc, argv, "!"); + + if (inputFile == "!") { + cout << "Missing file specification. Aborting. (-h for help)" << endl; + return -1; + } + + ifstream inFile(inputFile.c_str()); + if (!inFile) { + ERROR("File not found.", exception); + } + + CifLoader il(inFile); + il.setNoHAtoms(); + allCh = il.getAllChains(); + + for (unsigned int i = 0; i < allCh.size(); i++) { + cout << "\t," << allCh[i] << ","; + } + cout << "\n"; + + /*check on validity of chain: + if user select a chain then check validity + else select first valid one by default*/ + if (chainID != "!") { + bool validChain = false; + for (unsigned int i = 0; i < allCh.size(); i++) { + if (allCh[i] == chainID[0]) { + il.setChain(chainID[0]); + cout << "Loading chain " << chainID << "\n"; + validChain = true; + break; + } + } + if (!validChain) { + cout << "Chain " << chainID << " is not available\n"; + return -1; + } + + } else { + chainID[0] = allCh[0]; + cout << "Using chain " << chainID << "\n"; + } + + Protein prot; + prot.load(il); + Spacer *sp; + sp = prot.getSpacer(chainID[0]); + + allCh = il.getAllChains(); + cout << ">" << inputFile << "\n"; + + inFile.close(); + + for (unsigned int i = 0; i < sp->sizeAmino(); i++) { + switch (sp->getAmino(i).getState()) { + case HELIX: + cout << "H"; + break; + case STRAND: + cout << "E"; + break; + default: + cout << "."; + }; + if ((i + 1) % 60 == 0) + cout << "\n"; + } + cout << "\n"; + + return 0; +} + diff --git a/Biopool/APPS/Cif2Seq.cc b/Biopool/APPS/Cif2Seq.cc new file mode 100644 index 0000000..b62e096 --- /dev/null +++ b/Biopool/APPS/Cif2Seq.cc @@ -0,0 +1,133 @@ +/* This file is part of Victor. + + Victor is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + Victor 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 General Public License for more details. + + You should have received a copy of the GNU General Public License + along with Victor. If not, see . + */ + +#include +#include +#include +#include +#include + +using namespace Victor; +using namespace Victor::Biopool; + +void sShowHelp() { + cout << "Cif 2 Seq $Revision: 0.1 $ -- converts a CIF file into SEQ\n" + << "(torsion angles) protein structure backbone torsion angles\n" + << " Options: \n" + << "\t-i \t Input CIF file\n" + << "\t-o \t Output to file (default stdout)\n" + << "\t-c \t Chain identifier to read\n" + << "\t--all \t All chains\n" + << "\t-m \t Model number to read (NMR only, default is first model)\n" + << "\t--chi \t Write Chi angles (default false)\n" + << "\t-v \t verbose output\n\n" + << "\tIf both -c and --all are missing, only the first chain is processed.\n\n"; + +} + +int main(int argc, char** argv) { + + if (getArg("h", argc, argv)) { + sShowHelp(); + return 1; + } + + string inputFile, outputFile, chainID; + unsigned int modelNum; + bool chi, all; + + getArg("i", inputFile, argc, argv, "!"); + getArg("o", outputFile, argc, argv, "!"); + getArg("c", chainID, argc, argv, "!"); + getArg("m", modelNum, argc, argv, 999); + all = getArg("-all", argc, argv); + chi = getArg("-chi", argc, argv); + + // Check input file + if (inputFile == "!") { + cout << "Missing input file specification. Aborting. (-h for help)" << endl; + return -1; + } + ifstream inFile(inputFile.c_str()); + if (!inFile) + ERROR("Input file not found.", exception); + + + CifLoader cl(inFile); + + // Set PdbLoader variables + cl.setModel(modelNum); + cl.setNoHAtoms(); + cl.setNoHetAtoms(); + cl.setNoSecondary(); + if (!getArg("v", argc, argv)) { + cl.setNoVerbose(); + } + + + // Check chain args + if ((chainID != "!") && all) { + ERROR("You can use --all or -c, not both", error); + } + // User selected chain + if (chainID != "!") { + if (chainID.size() > 1) + ERROR("You can choose only 1 chain", error); + cl.setChain(chainID[0]); + }// All chains + else if (all) { + cl.setAllChains(); + }// First chain + else { + cl.setChain(cl.getAllChains()[0]); + } + + // Load the protein object + Protein prot; + prot.load(cl); + + inFile.close(); + + // Open the proper output stream (file or stdout) + std::ostream* os = &cout; + std::ofstream fout; + if (outputFile != "!") { + fout.open(outputFile.c_str()); + if (!fout) { + ERROR("Could not open file for writing.", exception); + } else { + os = &fout; + } + } + + + Spacer* sp; + for (unsigned int i = 0; i < prot.sizeProtein(); i++) { + + sp = prot.getSpacer(i); + + // Write the sequence + SeqSaver ss(*os); + if (!chi) + ss.setWriteChi(false); + sp->save(ss); + } + + fout.close(); + + return 0; +} + diff --git a/Biopool/APPS/CifCorrector.cc b/Biopool/APPS/CifCorrector.cc new file mode 100644 index 0000000..0882727 --- /dev/null +++ b/Biopool/APPS/CifCorrector.cc @@ -0,0 +1,71 @@ +/* This file is part of Victor. + + Victor is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + Victor 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 General Public License for more details. + + You should have received a copy of the GNU General Public License + along with Victor. If not, see . + */ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace Victor; +using namespace Victor::Biopool; + +/* + * + */ +int main(int argc, char** argv) { + + if (argc != 2) { + cout << "Cif Corrector $Revision: 0.1 $ -- adds missing oxygen atoms to " + << "protein structure backbones" << endl; + cout << " Usage: \t\t CifCorrector \n"; + return 1; + }; + + ifstream inFile(argv[1]); + if (!inFile) + ERROR("File not found.", exception); + + CifLoader il(inFile); + Protein prot; + prot.load(il); + unsigned int zero = 0; + Spacer sp = *(prot.getSpacer(zero)); + + for (unsigned int i = 0; i < sp.sizeAmino(); i++) + sp.getAmino(i).addMissingO(); + + inFile.close(); + + ofstream outFile(argv[1]); + + if (!outFile) + ERROR("Couldn't write file.", exception); + + CifSaver pss(outFile); + sp.save(pss); + + outFile.close(); + + return 0; +} + diff --git a/Biopool/APPS/CifEditor.cc b/Biopool/APPS/CifEditor.cc new file mode 100644 index 0000000..39f7afd --- /dev/null +++ b/Biopool/APPS/CifEditor.cc @@ -0,0 +1,86 @@ +/* This file is part of Victor. + + Victor is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + Victor 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 General Public License for more details. + + You should have received a copy of the GNU General Public License + along with Victor. If not, see . + */ + +#include +#include + +using namespace Victor; +using namespace Victor::Biopool; + +int main(int argc, char** argv) { + if (argc != 3) { + cout << "CIF Editor $Revision: 0.1 $ -- allows sequential manipulation of " + << "protein structure backbone torsion angles" << endl; + cout << " Usage: \t\t CifEditor \n"; + return 1; + } + + ifstream inFile(argv[1]); + if (!inFile) + ERROR("Input file not found.", exception); + + CifLoader cl(inFile); + Protein prot; + prot.load(cl); + unsigned int zero = 0; + Spacer sp = *(prot.getSpacer(zero)); + + cout << "Editing " << argv[1] << " output goes to " << argv[2] << "\n"; + + int aaid = -1; + do { + cout << "Aminoacid# or -1: "; + cin >> aaid; + if (aaid <= -1) { + cout << "Bye.\n"; + return 0; + } + + if (aaid >= (int) sp.sizeAmino()) { + cout << "\t Invalid aa#!\n"; + } else { + double newVal = 999; + cout << "\t " << aaid << " " << sp.getAmino(aaid).getType() << "\n"; + cout << "Phi= " << sp.getAmino(aaid).getPhi() + << "\t new phi or 999: "; + cin >> newVal; + if (newVal != 999) + sp.getAmino(aaid).setPhi(newVal); + cout << "Psi= " << sp.getAmino(aaid).getPsi() + << "\t new psi or 999: "; + cin >> newVal; + if (newVal != 999) + sp.getAmino(aaid).setPsi(newVal); + cout << "Omega= " << sp.getAmino(aaid).getOmega() + << "\t new omega or 999: "; + cin >> newVal; + if (newVal != 999) + sp.getAmino(aaid).setOmega(newVal); + + ofstream outFile2(argv[2]); + + if (!outFile2) + ERROR("Couldn't write output file.", exception); + + CifSaver pss2(outFile2); + + sp.save(pss2); + } + } while (aaid != -1); + + return 0; +} + diff --git a/Biopool/APPS/CifMover.cc b/Biopool/APPS/CifMover.cc new file mode 100644 index 0000000..74dc560 --- /dev/null +++ b/Biopool/APPS/CifMover.cc @@ -0,0 +1,204 @@ +/* This file is part of Victor. + + Victor is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + Victor 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 General Public License for more details. + + You should have received a copy of the GNU General Public License + along with Victor. If not, see . + */ + +#include +#include +#include +#include +#include +#include +#include + +using namespace Victor; +using namespace Victor::Biopool; + +// minimum distance between neighbouring CAs +const double LAMBDA = 1.5; + +void sShowHelp() { + cout << "CIF Mover\n" + << "Allows to move all residues in file by fixed offset.\n" + << " Options: \n" + << "\t-i \t\t Input CIF file\n" + << "\t-o \t\t Output CIF file\n" + << "\t[-r ] \t\t Repeat unit length (default = no rotation)\n" + << "\t[-l ] \t\t Angle lambda factor (default = 0.1)\n" + << "\t[-s ] \t\t Start residue of fragment (default = first)\n" + << "\t[-e ] \t\t End residue of fragment (default = last)\n" + << "\n"; +} + +void sAddLine() { + cout << "-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-\n"; +} + +int main(int argc, char** argv) { + + // -------------------------------------------------- + // 0. treat options + // -------------------------------------------------- + + if (getArg("h", argc, argv)) { + sShowHelp(); + return 1; + }; + + string inputFile, outputFile; + unsigned int startOffset, endOffset, repeatLength; + double lambdaAngle; + + getArg("i", inputFile, argc, argv, "!"); + getArg("o", outputFile, argc, argv, "!"); + getArg("s", startOffset, argc, argv, 0); + getArg("e", endOffset, argc, argv, 9999); + getArg("r", repeatLength, argc, argv, 9999); + getArg("l", lambdaAngle, argc, argv, 0.1); + + vgVector3 transOff; + for (unsigned int i = 0; i < 3; i++) + transOff[i] = 0.0; + + if ((inputFile == "!") || (outputFile == "!")) { + cout << "Missing file specification. Aborting. (-h for help)" << endl; + return -1; + } + + // -------------------------------------------------- + // 1. read structure + // -------------------------------------------------- + + ifstream inFile(inputFile.c_str()); + + if (!inFile) + ERROR("File does not exist.\n", exception); + + CifLoader cl(inFile); + Protein prot; + prot.load(cl); + unsigned int zero = 0; + Spacer sp = *(prot.getSpacer(zero)); + + inFile.close(); + + + endOffset = sp.getIndexFromPdbNumber(endOffset); + if (startOffset > 0) + startOffset = sp.getIndexFromPdbNumber(startOffset); + + // -------------------------------------------------- + // 2. rotate spacer + // -------------------------------------------------- + + if (repeatLength < 9999) { + IntCoordConverter icc; + + // find offset + vgVector3 firstA = sp.getAmino(endOffset + - (3 * repeatLength / 4))[CA].getCoords() + - sp.getAmino(endOffset)[CA].getCoords(); + + vgVector3 firstB = sp.getAmino(endOffset + - (1 * repeatLength / 4))[CA].getCoords() + - sp.getAmino(endOffset)[CA].getCoords(); + + vgVector3 firstNorm = (firstA.normalize()).cross(firstB.normalize()); + + vgVector3 secondA = sp.getAmino(endOffset + repeatLength + - (3 * repeatLength / 4))[CA].getCoords() + - sp.getAmino(endOffset)[CA].getCoords(); + + vgVector3 secondB = sp.getAmino(endOffset + repeatLength + - (1 * repeatLength / 4))[CA].getCoords() + - sp.getAmino(endOffset)[CA].getCoords(); + + vgVector3 secondNorm = (secondA.normalize()).cross(secondB.normalize()); + + firstNorm.normalize(); + secondNorm.normalize(); + + double scalar = icc.getAngle(firstNorm, secondNorm) * lambdaAngle; + + cout << "Scalar = " << setw(5) << setprecision(3) + << RAD2DEG * scalar / lambdaAngle << "\n"; + + + vgVector3 axis = (firstNorm.normalize()).cross(secondNorm.normalize()); + vgMatrix3 res = vgMatrix3::createRotationMatrix(axis, scalar); + + sp.getAmino(startOffset)[N].addRot(res); + } + + // -------------------------------------------------- + // 3. translate spacer + // -------------------------------------------------- + + if (endOffset < sp.sizeAmino() - 1) { + sp.getAmino(endOffset).unbindOut(sp.getAmino(endOffset + 1)); + sp.getAmino(endOffset)[C].unbindOut(sp.getAmino(endOffset + 1)[N]); + + // find offset + vgVector3 first = sp.getAmino(endOffset)[C].getCoords(); + vgVector3 second = sp.getAmino(endOffset + 1)[N].getCoords(); + + double d = sp.getAmino(endOffset)[C].distance( + sp.getAmino(endOffset + 1)[N]); + + double frac = (d - LAMBDA) / d; + + for (unsigned int i = 0; i < 3; i++) + transOff[i] = (second[i] - first[i]) * frac; + } else { + if (startOffset != 0) + ERROR("Both start and end offset are undefined.", exception); + + endOffset = sp.sizeAmino() - 1; + + // NB: unbind has to be reversed after moving the atoms if the model + // is to be used further in the same program + sp.getAmino(startOffset - 1).unbindOut(sp.getAmino(startOffset)); + sp.getAmino(startOffset - 1)[C].unbindOut(sp.getAmino(startOffset)[N]); + + // find offset + vgVector3 first = sp.getAmino(startOffset + 1)[N].getCoords(); + vgVector3 second = sp.getAmino(startOffset)[C].getCoords(); + + double d = sp.getAmino(startOffset + 1)[N].distance( + sp.getAmino(startOffset)[C]); + + double frac = (d - LAMBDA) / d; + + for (unsigned int i = 0; i < 3; i++) + transOff[i] = (second[i] - first[i]) * frac; + } + + sp.getAmino(startOffset)[N].addTrans(transOff); + + // -------------------------------------------------- + // 4. write model to disk + // -------------------------------------------------- + + ofstream outFile(outputFile.c_str()); + if (!outFile) { + ERROR("File not found.", exception); + } + CifSaver cs(outFile); + sp.save(cs); + + outFile.close(); + + return 0; +} + diff --git a/Biopool/APPS/CifReaderWriter.cc b/Biopool/APPS/CifReaderWriter.cc new file mode 100644 index 0000000..85487d6 --- /dev/null +++ b/Biopool/APPS/CifReaderWriter.cc @@ -0,0 +1,120 @@ +/* This file is part of Victor. + + Victor is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + Victor 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 General Public License for more details. + + You should have received a copy of the GNU General Public License + along with Victor. If not, see . + */ + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace Victor; +using namespace Victor::Biopool; +using namespace std; + +void sShowHelp() { + cout << "CIF Reader / Writer \n" + << "Allows to read and write a CIF or a PDB file without modifying them.\n" + << " Options: \n" + << "\t-i \t\t Input CIF/PDB file\n" + << "\t-o \t\t Output CIF/PDB file\n" + << "\n"; +} + +int main(int argc, char** argv) { + + // -------------------------------------------------- + // 0. treat options + // -------------------------------------------------- + + if (getArg("h", argc, argv)) { + sShowHelp(); + return 1; + }; + + string inputFile, outputFile; + Protein prot; + + getArg("i", inputFile, argc, argv, "!"); + getArg("o", outputFile, argc, argv, "!"); + + if ((inputFile == "!") || (outputFile == "!")) { + cout << "Missing file specification. Aborting. (-h for help)" << endl; + return -1; + } + + // -------------------------------------------------- + // 1. read structure + // -------------------------------------------------- + + ifstream inFile(inputFile.c_str()); + if (!inFile) { + ERROR("File does not exist.\n", exception); + } + + if (inputFile.find("pdb") != string::npos) { + cout << "Loading PDB file..." << endl; + + PdbLoader pl(inFile); + prot.load(pl); + + } else if (inputFile.find("cif") != string::npos) { + cout << "Loading CIF file..." << endl; + + CifLoader cl(inFile); + prot.load(cl); + + } else { + cout << "Uknown input file format. Aborting. (-h for help)" << endl; + return -2; + } + + inFile.close(); + + // -------------------------------------------------- + // 2. write structure + // -------------------------------------------------- + + ofstream outFile(outputFile.c_str()); + if (!outFile) { + ERROR("File not found.", exception); + } + + if (outputFile.find("pdb") != string::npos) { + cout << "Saving PDB file..." << endl; + + PdbSaver ps(outFile); + prot.save(ps); + + } else if (outputFile.find("cif") != string::npos) { + cout << "Saving CIF file..." << endl; + + CifSaver cs(outFile); + prot.save(cs); + + } else { + cout << "Uknown output file format. Aborting. (-h for help)" << endl; + return -3; + } + + outFile.close(); + + return 0; +} + diff --git a/Biopool/APPS/CifSecondary.cc b/Biopool/APPS/CifSecondary.cc new file mode 100644 index 0000000..cde7b87 --- /dev/null +++ b/Biopool/APPS/CifSecondary.cc @@ -0,0 +1,221 @@ +/* This file is part of Victor. + + Victor is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + Victor 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 General Public License for more details. + + You should have received a copy of the GNU General Public License + along with Victor. If not, see . + */ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace Victor; +using namespace Victor::Biopool; + +void sShowHelp() { + cout << "CIF Secondary $Revision: 0.1 $ -- calculate the secondary structure\n" + << "(torsion angles) protein structure backbone torsion angles\n" + << "\tOptions: \n" + << "\t-i \t Input CIF file\n" + << "\t-o \t Output to file(the chain letter is appended)\n" + << "\t-c \t Chain identifier to read(default is all chains)\n" + << "\t-m \t Model number to read (NMR only, default is first model)\n" + << "\t-s <1,2,3> \t SS calculation(default 3): 1 = CIF fields, 2 = torsion angles, 3 = DSSP\n" + << "\t \t 1,2:\t H = helix, E = extended(strand,sheet), . = other.\n" + << "\t \t 3:\t H = alpha-helix, E = sheet, B = bridge,\n" + << "\t \t G = 3-10-helix, I = pi-helix, T = n-Turn, S = bend\n" + << "\t--ext \t Extended output (default false). Write the type of aminoacid:\n" + << "\t \t N = negative charge, \t P = positive charge\n" + << "\t \t h = hydrophilic, \t + = hydrophobic\n" + << "\t \t , = neutral (hydrophobic)\n" + << "\t-v \t verbose output\n\n" + << "\tWhen parsing multiple chains, one file for each chain is created\n\n"; +} + +void writeOutput(Spacer* sp, bool ext, int ssType, ostream& os) { + + if (ext) { + for (unsigned int i = 0; i <= sp->sizeAmino(); i++) { + if ((i + 1) % 10 == 0) + os << setw(1) << (((i + 1) % 100) / 10); + else + os << " "; + if ((i + 1) % 60 == 0) + os << "\n"; + } + os << "\n"; + + for (unsigned int i = 0; i < sp->sizeAmino(); i++) { + os << sp->getAmino(i).getType1L(); + if ((i + 1) % 60 == 0) + os << "\n"; + } + os << "\n"; + + for (unsigned int i = 0; i < sp->sizeAmino(); i++) { + switch (sp->getAmino(i).getCode()) { + case ASP: + case GLU: + os << "P"; + break; + case LYS: + case ARG: + os << "N"; + break; + case ASN: + case GLN: + case SER: + case THR: + case HIS: + os << "h"; + break; + case VAL: + case LEU: + case ILE: + os << "+"; + break; + default: + os << ","; + }; + if ((i + 1) % 60 == 0) + os << "\n"; + } + os << "\n"; + } + if (ssType != 3) { + for (unsigned int i = 0; i < sp->sizeAmino(); i++) { + switch (sp->getAmino(i).getState()) { + case HELIX: + os << "H"; + break; + case STRAND: + os << "E"; + break; + default: + os << "."; + }; + if ((i + 1) % 60 == 0) + os << "\n"; + } + } else { + vector > ss = sp->getDSSP(); + for (unsigned int i = 0; i < ss.size(); i++) { + if (!ss[i].empty()) { + for (set ::iterator it = ss[i].begin(); it != ss[i].end(); ++it) { + os << (*it); + } + } else { + os << "."; + } + if ((i + 1) % 60 == 0) + os << "\n"; + } + } + os << "\n"; +} + +int main(int argc, char** argv) { + + if (getArg("h", argc, argv)) { + sShowHelp(); + return 1; + } + + string inputFile, outputFile, chainID; + unsigned int modelNum; + unsigned int ssType; + bool extendedOutput, all; + + getArg("i", inputFile, argc, argv, "!"); + getArg("o", outputFile, argc, argv, "!"); + getArg("c", chainID, argc, argv, "!"); + getArg("m", modelNum, argc, argv, 999); + getArg("s", ssType, argc, argv, 3); + all = getArg("-all", argc, argv); + extendedOutput = getArg("-ext", argc, argv); + + // Check input file + if (inputFile == "!") { + cout << "Missing input file specification. Aborting. (-h for help)" << endl; + return -1; + } + ifstream inFile(inputFile.c_str()); + if (!inFile) + ERROR("Input file not found.", exception); + + CifLoader cl(inFile); + + // Set CifLoader variables + cl.setModel(modelNum); + cl.setNoHetAtoms(); + + if (!getArg("v", argc, argv)) { + cl.setNoVerbose(); + } + + // Check chain args + if ((chainID != "!") && all) { + ERROR("You can use --all or -c, not both", error); + } + // User selected chain + if (chainID != "!") { + if (chainID.size() > 1) + ERROR("You can choose only 1 chain", error); + cl.setChain(chainID[0]); + }// All chains + else if (all) { + cl.setAllChains(); + }// First chain + else { + cl.setChain(cl.getAllChains()[0]); + } + + // Load the protein object + Protein prot; + prot.load(cl); + + inFile.close(); + + // Open the proper output stream (file or stdout) + std::ostream* os = &cout; + std::ofstream fout; + if (outputFile != "!") { + fout.open(outputFile.c_str()); + if (!fout) { + ERROR("Could not open file for writing.", exception); + } else { + os = &fout; + } + } + + Spacer* sp; + for (unsigned int i = 0; i < prot.sizeProtein(); i++) { + + sp = prot.getSpacer(i); + writeOutput(sp, extendedOutput, ssType, (*os)); + } + + fout.close(); + + return 0; +} + diff --git a/Biopool/APPS/CifShifter.cc b/Biopool/APPS/CifShifter.cc new file mode 100644 index 0000000..e620c64 --- /dev/null +++ b/Biopool/APPS/CifShifter.cc @@ -0,0 +1,142 @@ +/* This file is part of Victor. + + Victor is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + Victor 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 General Public License for more details. + + You should have received a copy of the GNU General Public License + along with Victor. If not, see . + */ + +#include +#include +#include +#include + +using namespace Victor; +using namespace Victor::Biopool; + +void sShowHelp() { + cout << "CIF Shifter\n" + << "Allows to shift all residues in file by fixed offset.\n" + << " Options: \n" + << "\t-i \t\t Input CIF file\n" + << "\t-o \t\t Output CIF file\n" + << "\t[-p ] \t\t Positive *residue* offset\n" + << "\t[-n ] \t\t Negative *residue* offset\n" + << "\t[-P ] \t\t Positive *atom* offset\n" + << "\t[-N ] \t\t Negative *atom* offset\n" + << "\t[--nohydrogen] \t\t Skip hydrogen atoms\n" + << "\t[--renum] \t\t Reset residue numbering starting from 1\n" + << "\n"; +} + +void sAddLine() { + cout << "-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-\n"; +} + +void sRenumberAtoms(Spacer& sp) { + unsigned int counter = 1; + for (unsigned int i = 0; i < sp.sizeAmino(); i++) + for (unsigned int j = 0; j < sp.getAmino(i).size(); j++) { + sp.getAmino(i)[j].setNumber(counter); + counter++; + } +} + +int main(int argc, char** argv) { + + // -------------------------------------------------- + // 0. treat options + // -------------------------------------------------- + + if (getArg("h", argc, argv)) { + sShowHelp(); + return 1; + }; + + string inputFile, outputFile; + int offset, offsetAtom, tmp; + getArg("i", inputFile, argc, argv, "!"); + getArg("o", outputFile, argc, argv, "!"); + + getArg("p", offset, argc, argv, 0); + getArg("n", tmp, argc, argv, 0); + offset -= tmp; + + getArg("P", offsetAtom, argc, argv, 0); + getArg("N", tmp, argc, argv, 0); + offsetAtom -= tmp; + + bool noHydrogen = getArg("-nohydrogen", argc, argv); + bool renumber = getArg("-renum", argc, argv); + + if ((inputFile == "!") || (outputFile == "!")) { + cout << "Missing file specification. Aborting. (-h for help)" << endl; + return -1; + } + + if ((offset == 0) && (!renumber) && (!noHydrogen)) { + cout << "Warning: Offset is zero. Mistake? \n"; + } + // -------------------------------------------------- + // 1. read structure + // -------------------------------------------------- + + ifstream inFile(inputFile.c_str()); + + if (!inFile) + ERROR("File does not exist.\n", exception); + + CifLoader cl(inFile); + + if (noHydrogen) + cl.setNoHAtoms(); + + Protein prot; + prot.load(cl); + unsigned int zero = 0; + Spacer sp = *(prot.getSpacer(zero)); + inFile.close(); + + // -------------------------------------------------- + // 1.1 renumber residues (if necessary) + // -------------------------------------------------- + + if (renumber) { + cout << "Renumbering...\n"; + sp.setStartOffset(1); + sp.removeAllGaps(); + sRenumberAtoms(sp); + } + + // -------------------------------------------------- + // 2. shift offset + // -------------------------------------------------- + + sp.setStartOffset(offset + sp.getStartOffset()); + + if (offsetAtom != 0) + sp.setAtomStartOffset(offsetAtom + sp.getAtomStartOffset()); + + // -------------------------------------------------- + // 3. write model to disk + // -------------------------------------------------- + + ofstream outFile(outputFile.c_str()); + if (!outFile) + ERROR("File not found.", exception); + CifSaver cs(outFile); + + sp.save(cs); + outFile.close(); + + return 0; +} + diff --git a/Biopool/APPS/Makefile b/Biopool/APPS/Makefile index cb26b82..d52c2b6 100644 --- a/Biopool/APPS/Makefile +++ b/Biopool/APPS/Makefile @@ -27,16 +27,20 @@ INC_PATH = -I. -I../../tools/ -I../../Biopool/Sources -I../../Energy/Sources -I. # SOURCES = PdbCorrector.cc PdbSecondary.cc PdbEditor.cc Pdb2Seq.cc pdb2secondary.cc pdbshifter.cc \ - pdbMover.cc + pdbMover.cc CifEditor.cc CifSecondary.cc Cif2Secondary.cc CifMover.cc \ + CifShifter.cc CifCorrector.cc Cif2Seq.cc CifReaderWriter.cc OBJECTS = PdbCorrector.o PdbSecondary.o PdbEditor.o Pdb2Seq.o pdb2secondary.o pdbshifter.o \ - pdbMover.o + pdbMover.o CifEditor.o CifSecondary.o Cif2Secondary.o CifMover.o \ + CifShifter.o CifCorrector.o Cif2Seq.o CifReaderWriter.o TARGETS = PdbCorrector PdbSecondary PdbEditor Pdb2Seq pdb2secondary pdbshifter \ - pdbMover + pdbMover CifEditor CifSecondary Cif2Secondary CifMover CifShifter \ + CifCorrector Cif2Seq CifReaderWriter EXECS = PdbCorrector PdbSecondary PdbEditor Pdb2Seq pdb2secondary pdbshifter \ - pdbMover + pdbMover CifEditor CifSecondary Cif2Secondary CifMover CifShifter \ + CifCorrector Cif2Seq CifReaderWriter LIBRARY = APPSlibBiopool.a diff --git a/Biopool/APPS/PdbCorrector.cc b/Biopool/APPS/PdbCorrector.cc index 83e6b04..777245f 100644 --- a/Biopool/APPS/PdbCorrector.cc +++ b/Biopool/APPS/PdbCorrector.cc @@ -36,16 +36,20 @@ int main(int nArgs, char* argv[]) { return 1; }; - Spacer sp; ifstream inFile(argv[1]); if (!inFile) ERROR("File not found.", exception); PdbLoader il(inFile); - sp.load(il); + Protein prot; + prot.load(il); + unsigned int zero = 0; + Spacer sp = *(prot.getSpacer(zero)); for (unsigned int i = 0; i < sp.sizeAmino(); i++) sp.getAmino(i).addMissingO(); + + inFile.close(); ofstream outFile2(argv[1]); @@ -55,5 +59,7 @@ int main(int nArgs, char* argv[]) { PdbSaver pss2(outFile2); sp.save(pss2); + outFile2.close(); + return 0; } diff --git a/Biopool/APPS/PdbEditor.cc b/Biopool/APPS/PdbEditor.cc index 33ab8fd..a952ee4 100644 --- a/Biopool/APPS/PdbEditor.cc +++ b/Biopool/APPS/PdbEditor.cc @@ -31,13 +31,15 @@ int main(int nArgs, char* argv[]) { return 1; }; - Spacer sp; ifstream inFile(argv[1]); if (!inFile) ERROR("File not found.", exception); PdbLoader il(inFile); - sp.load(il); + Protein prot; + prot.load(il); + unsigned int zero = 0; + Spacer sp = *(prot.getSpacer(zero)); cout << "Editing " << argv[1] << " output goes to " << argv[2] << "\n"; diff --git a/Biopool/APPS/pdbMover.cc b/Biopool/APPS/pdbMover.cc index 21d6b35..95a62f2 100644 --- a/Biopool/APPS/pdbMover.cc +++ b/Biopool/APPS/pdbMover.cc @@ -78,16 +78,16 @@ int main(int nArgs, char* argv[]) { // 1. read structure // -------------------------------------------------- - Spacer sp; - ifstream inFile(inputFile.c_str()); if (!inFile) ERROR("File does not exist.\n", exception); PdbLoader pl(inFile); - - sp.load(pl); + Protein prot; + prot.load(pl); + unsigned int zero = 0; + Spacer sp = *(prot.getSpacer(zero)); inFile.close(); diff --git a/Biopool/APPS/pdbshifter.cc b/Biopool/APPS/pdbshifter.cc index 8a6d3db..07891d3 100644 --- a/Biopool/APPS/pdbshifter.cc +++ b/Biopool/APPS/pdbshifter.cc @@ -87,8 +87,6 @@ int main(int nArgs, char* argv[]) { // 1. read structure // -------------------------------------------------- - Spacer sp; - ifstream inFile(inputFile.c_str()); if (!inFile) @@ -99,7 +97,10 @@ int main(int nArgs, char* argv[]) { if (noHydrogen) pl.setNoHAtoms(); - sp.load(pl); + Protein prot; + prot.load(pl); + unsigned int zero = 0; + Spacer sp = *(prot.getSpacer(zero)); inFile.close(); // -------------------------------------------------- diff --git a/Biopool/Sources/AminoAcid.h b/Biopool/Sources/AminoAcid.h index f497833..0fd9e0e 100644 --- a/Biopool/Sources/AminoAcid.h +++ b/Biopool/Sources/AminoAcid.h @@ -58,9 +58,9 @@ namespace Victor { namespace Biopool { unsigned int size() const; unsigned int sizeBackbone() const; - double getPhi(bool override = false); - double getPsi(bool override = false); - double getOmega(bool override = false); + double getPhi(bool bypass = false); + double getPsi(bool bypass = false); + double getOmega(bool bypass = false); double getChi(unsigned int n); vector getChi(); unsigned int getMaxChi(); @@ -187,8 +187,8 @@ namespace Victor { namespace Biopool { } inline double - AminoAcid::getPhi(bool override) { - if ((phi > 990) || (override)) + AminoAcid::getPhi(bool bypass) { + if ((phi > 990) || (bypass)) if (sizeInBonds()) phi = RAD2DEG * icc.getTorsionAngle(getInBond(0)[C], (*this)[N], (*this)[CA], (*this)[C]); @@ -196,8 +196,8 @@ namespace Victor { namespace Biopool { } inline double - AminoAcid::getPsi(bool override) { - if ((psi > 990) || (override)) + AminoAcid::getPsi(bool bypass) { + if ((psi > 990) || (bypass)) if (sizeOutBonds()) psi = RAD2DEG * icc.getTorsionAngle((*this)[N], (*this)[CA], (*this)[C], getOutBond(0)[N]); @@ -205,8 +205,8 @@ namespace Victor { namespace Biopool { } inline double - AminoAcid::getOmega(bool override) { - if ((omega > 990) || (override)) + AminoAcid::getOmega(bool bypass) { + if ((omega > 990) || (bypass)) if (sizeOutBonds()) omega = RAD2DEG * icc.getTorsionAngle((*this)[CA], (*this)[C], getOutBond(0)[N], getOutBond(0)[CA]); diff --git a/Biopool/Sources/AminoAcidHydrogen.cc b/Biopool/Sources/AminoAcidHydrogen.cc index 9fd5003..4b6cc4f 100644 --- a/Biopool/Sources/AminoAcidHydrogen.cc +++ b/Biopool/Sources/AminoAcidHydrogen.cc @@ -121,6 +121,11 @@ AminoAcidHydrogen::setHydrogen(AminoAcid* aa, bool verbose) { IntCoordConverter icc; Atom atH; atH.setType("H"); + + // extra CIF fields + atH.setAsymId(aa->getAtom(0).getAsymId()); + atH.setEntityId(aa->getAtom(0).getEntityId()); + atH.setModel(aa->getAtom(0).getModel()); AminoAcid before = (*aa).getInBond(0); atH.bindIn((*aa)[N]); @@ -170,6 +175,11 @@ AminoAcidHydrogen::setHydrogen(AminoAcid* aa, bool verbose) { chiral = atoi(args[8].c_str()); atH.setType(args[2]); + + // extra CIF fields + atH.setAsymId(aa->getAtom(0).getAsymId()); + atH.setEntityId(aa->getAtom(0).getEntityId()); + atH.setModel(aa->getAtom(0).getModel()); if (aa->getSideChain().isMember(atBindCod)) { diff --git a/Biopool/Sources/Atom.cc b/Biopool/Sources/Atom.cc index ebc2f10..48fffc5 100644 --- a/Biopool/Sources/Atom.cc +++ b/Biopool/Sources/Atom.cc @@ -34,7 +34,7 @@ using namespace Victor; using namespace Victor::Biopool; */ Atom::Atom(unsigned int mI, unsigned int mO) : SimpleBond(mI, mO), superior(NULL), type(X), coords(0, 0, 0), Bfac(0.0), trans(0, 0, 0), rot(1), -modified(false) { +modified(false), asymId('X'), entityId("0"), occupancy(0.0), model(0) { PRINT_NAME; } @@ -133,6 +133,11 @@ Atom::copy(const Atom& orig) { type = orig.type; coords = orig.coords; Bfac = orig.Bfac; + + asymId = orig.asymId; + entityId = orig.entityId; + occupancy = orig.occupancy; + model = orig.model; trans = orig.trans; rot = orig.rot; diff --git a/Biopool/Sources/Atom.h b/Biopool/Sources/Atom.h index 3df6f47..2c1fa6c 100644 --- a/Biopool/Sources/Atom.h +++ b/Biopool/Sources/Atom.h @@ -53,6 +53,11 @@ namespace Victor { namespace Biopool { double getBFac() { return Bfac; } + + char getAsymId(); + string getEntityId(); + double getOccupancy(); + int getModel(); double distance(Atom& other); @@ -83,6 +88,10 @@ namespace Victor { namespace Biopool { void setBFac(double _b) { Bfac = _b; } + void setAsymId(char aId); + void setEntityId(string eId); + void setOccupancy(double occ); + void setModel(int mod); void setTrans(vgVector3 t); void addTrans(vgVector3 t); @@ -117,6 +126,11 @@ namespace Victor { namespace Biopool { vgVector3 coords; // xyz-Coords double Bfac; // B-factor + + char asymId; + string entityId; + double occupancy; + int model; vgVector3 trans; // relative translation vgMatrix3 rot; // relative rotation @@ -179,6 +193,22 @@ namespace Victor { namespace Biopool { Atom::inSync() { return (!modified); } + + inline char Atom::getAsymId() { + return asymId; + } + + inline string Atom::getEntityId() { + return entityId; + } + + inline double Atom::getOccupancy() { + return occupancy; + } + + inline int Atom::getModel() { + return model; + } // MODIFIERS: @@ -279,6 +309,22 @@ namespace Victor { namespace Biopool { this->superior = gr; } + inline void Atom::setAsymId(char aId) { + asymId = aId; + } + + inline void Atom::setEntityId(string eId) { + entityId = eId; + } + + inline void Atom::setOccupancy(double occ) { + occupancy = occ; + } + + inline void Atom::setModel(int mod) { + model = mod; + } + // OPERATORS: diff --git a/Biopool/Sources/CifLoader.cc b/Biopool/Sources/CifLoader.cc new file mode 100644 index 0000000..d4ad53a --- /dev/null +++ b/Biopool/Sources/CifLoader.cc @@ -0,0 +1,703 @@ +/* This file is part of Victor. + + Victor is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + Victor 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 General Public License for more details. + + You should have received a copy of the GNU General Public License + along with Victor. If not, see . + */ + +// Includes: +#include + +#include +#include +#include +#include +#include +#include +#include + +#include "CifLoader.h" + +// Global constants, typedefs, etc. (to avoid): + +using namespace Victor; +using namespace Victor::Biopool; +using namespace std; + +// CONSTRUCTORS/DESTRUCTOR: + +CifLoader::CifLoader(istream& _input, ostream& output, bool _permissive, bool _noHAtoms, + bool _noHetAtoms, bool _noSecondary, bool _noConnection, bool _noWater, + bool _verb, bool _allChains, string _NULL, bool _onlyMetal, + bool _noNucleotideChains ) : +input(_input), output(output), permissive(_permissive), valid(true), noHAtoms(_noHAtoms), +noHetAtoms(_noHetAtoms), noSecondary(_noSecondary), noConnection(_noConnection), +noWater(_noWater), verbose(_verb), allChains(_allChains), chain(' '), +model(999), altAtom('A'), helixCode(_NULL), +//sheetCode(_NULL), helixData(), sheetData(), onlyMetalHetAtoms(_onlyMetal), +sheetCode(_NULL), onlyMetalHetAtoms(_onlyMetal), noNucleotideChains(_noNucleotideChains) { + cif = new CifStructure(_input, output); +} + +CifLoader::~CifLoader() { + delete cif; + PRINT_NAME; +} + +// PREDICATES: + +/** + * If user selected a Model, it check validity of this choice, + * otherwise it select first available chain. + */ +void CifLoader::checkModel() { + if ((model != 999) && (model > getMaxModels())) { + ERROR("Please check model number", exception); + } +} + +/** + * If user selected a chain, it check validity of this choice, + * otherwise it select first available chain. + */ +void CifLoader::checkAndSetChain() { + vector chainList = getAllChains(); + + if (chain != ' ') { + bool validChain = false; + for (unsigned int i = 0; i < chainList.size(); i++) + if (chain == chainList[i]) { + validChain = true; + break; + } + if (validChain == false) { + ERROR("Please check chain id. This is not valid", exception); + } + } else { + chain = chainList[0]; //the first valid chain is default choice + } +} + +/** + * Reads in the maximum allowed number of NMR models, zero otherwise. + */ +unsigned int CifLoader::getMaxModels() { + input.clear(); // reset file error flags + input.seekg(0); + + string atomLine = readLine(input); + + unsigned int max = 0; + + // search column's number of the model field in the atom group + cif->parseGroup("atom", atomLine); + int col = cif->getGroupColumnNumber("atom", "model"); + + if (col != 0) { + while (input) { + if (atomLine.substr(0, 4) == "ATOM") { + max = stoiDEF(cif->getGroupField("atom", atomLine, col)); + } + atomLine = readLine(input); + } + } + return max; +} + +/** + * Returns all available chain IDs for a PDB file. + * @return all available chain IDs + */ +vector CifLoader::getAllChains() { + vector res; + char lastChain = ' '; + + input.clear(); // reset file error flags + input.seekg(0); + + string atomLine = readLine(input); + + unsigned int modelNum = 0; + + cif->parseGroup("atom", atomLine); + int modelCol = cif->getGroupColumnNumber("atom", "model"); + int chainCol = cif->getGroupColumnNumber("atom", "chain"); + + while (input) { + if (atomLine.substr(0, 4) == "ATOM") { + modelNum = stoiDEF(cif->getGroupField("atom", atomLine, modelCol));; + // only consider first model: others duplicate chain IDs + if (modelNum > 1) { + break; + } + // check for new chains containing amino acids + char id = (cif->getGroupField("atom", atomLine, chainCol).c_str())[0]; + if (id != lastChain) { + lastChain = id; + res.push_back(id); + } + } + atomLine = readLine(input); + } + return res; +} + +void CifLoader::setOnlyMetalHetAtoms() { + if (noHetAtoms) { + ERROR("can't load metal ions if hetAtoms option is disabled", exception); + } + onlyMetalHetAtoms = true; + noWater = true; +} + +void CifLoader::setWater() { + if (noHetAtoms || onlyMetalHetAtoms) { + ERROR("can't load water if hetAtoms option is disabled\nor onlyMetalHetAtoms is enabled", exception); + } + noWater = false; +} + +// HELPERS + +/** + * Private helper function to set bond structure after loading the spacer. + * @param Spacer reference + * @return bool + */ +bool CifLoader::setBonds(Spacer& sp) { + //cout << sp.getAmino(0).getType1L() << "\n"; + sp.getAmino(0).setBondsFromPdbCode(true); + for (unsigned int i = 1; i < sp.size(); i++) { + //cout << sp.getAmino(i).getType1L() << "\n"; + if (!sp.getAmino(i).setBondsFromPdbCode(true, &(sp.getAmino(i - 1)))) + return false; + } + return true; +} + +/** + * Private helper function to determine if atom is backbone or sidechain. + * @param Spacer reference + * @return bool + */ +bool CifLoader::inSideChain(const AminoAcid& aa, const Atom& at) { + if (isBackboneAtom(at.getCode())) { + return false; + } + if ((at.getType() == "H") || (at.getType() == "HN") + || ((at.getType() == "HA") && (!aa.isMember(HA))) + || (at.getType() == "1HA") || (at.getType() == "1H") + || (at.getType() == "2H") || (at.getType() == "3H")) { + return false; // special case for GLY H (code HA) + } + return true; // rest of aminoacid is its sidechain +} + +/** + * Try to assigns the secondary structure from the PDB header. If not present + * uses Spacer's setStateFromTorsionAngles(). + * @param Spacer reference + */ +void CifLoader::assignSecondary(Spacer& sp) { + if (helixData.size() + sheetData.size() == 0) { + sp.setStateFromTorsionAngles(); + return; + } + + for (unsigned int i = 0; i < helixData.size(); i++) { + if (helixCode[i] == chain) { + for (int j = helixData[i].first; j <= const_cast (helixData[i].second); j++) { + // important: keep ifs separated to avoid errors + if (j < sp.maxPdbNumber()) { + if (!sp.isGap(sp.getIndexFromPdbNumber(j))) { + sp.getAmino(sp.getIndexFromPdbNumber(j)).setState(HELIX); + } + } + } + } + } + + for (unsigned int i = 0; i < sheetData.size(); i++) { + if (sheetCode[i] == chain) { + for (int j = sheetData[i].first; j <= const_cast (sheetData[i].second); j++) { + // important: keep ifs separated to avoid errors + if (j < sp.maxPdbNumber()) { + if (!sp.isGap(sp.getIndexFromPdbNumber(j))) { + sp.getAmino(sp.getIndexFromPdbNumber(j)).setState(STRAND); + } + } + } + } + } +} + +/* +void +CifLoader::loadSpacer(Spacer& sp){ + Protein prot; + CifLoader::loadProtein(prot); + sp = prot.getSpacer(0); +} + */ + +/** + * Parse a single line of a CIF file. + * @param atomLine the whole CIF line as it is + * @param tag the first field (keyword) in a CIF line + * @param lig pointer to a ligan + * @param aa pointer to an amino acid + * @return Residue number read from the CIF line. + */ +int +CifLoader::parseCifline(string atomLine, string tag, Ligand* lig, AminoAcid* aa) { + // get atom id + int atNum = stoiDEF(cif->getGroupField("atom", atomLine, + cif->getGroupColumnNumber("atom", "atom id"))); + // get residue number + int aaNum = stoiDEF(cif->getGroupField("atom", atomLine, + cif->getGroupColumnNumber("atom", "residue num"))); + // get code for insertion of residues + char altAaID = cif->getGroupField("atom", atomLine, + cif->getGroupColumnNumber("atom", "residue ins")).c_str()[0]; + + // get x, y, z coordinates + vgVector3 coord; + coord.x = stodDEF(cif->getGroupField("atom", atomLine, + cif->getGroupColumnNumber("atom", "x"))); + coord.y = stodDEF(cif->getGroupField("atom", atomLine, + cif->getGroupColumnNumber("atom", "y"))); + coord.z = stodDEF(cif->getGroupField("atom", atomLine, + cif->getGroupColumnNumber("atom", "z"))); + + // get b-factor + double bfac = 0.0; + int colBfac = cif->getGroupColumnNumber("atom", "bfac"); + if (colBfac != -1) { + string sbfac = cif->getGroupField("atom", atomLine, colBfac); + if (sbfac != "?" && sbfac != ".") { + bfac = stodDEF(sbfac); + } + } + + // get atom name + string atType = cif->getGroupField("atom", atomLine, + cif->getGroupColumnNumber("atom", "atom name")); + + // get residue name + string aaType = cif->getGroupField("atom", atomLine, + cif->getGroupColumnNumber("atom", "residue name")); + + // take care of deuterium atoms + if (atType == "D") { + cerr << "--> " << atType << "\n"; + atType = "H"; + } + + // get asym id + char asymId = cif->getGroupField("atom", atomLine, + cif->getGroupColumnNumber("atom", "atom asym")).c_str()[0]; + + // get entity id + string entityId = cif->getGroupField("atom", atomLine, + cif->getGroupColumnNumber("atom", "atom entity")); + + // get occupancy + double occ = stodDEF(cif->getGroupField("atom", atomLine, + cif->getGroupColumnNumber("atom", "occupancy"))); + + // get model + int model = stoiDEF(cif->getGroupField("atom", atomLine, + cif->getGroupColumnNumber("atom", "model"))); + + // Initialize the Atom object + Atom* at = new Atom(); + at->setNumber(atNum); + at->setType(atType); + at->setCoords(coord); + at->setBFac(bfac); + at->setAsymId(asymId); + at->setEntityId(entityId); + at->setOccupancy(occ); + at->setModel(model); + + // Ligand object (includes DNA/RNA in "ATOM" field) + if (tag == "HETATM" || + isKnownNucleotide(nucleotideThreeLetterTranslator(aaType))) { + if (noWater) { + if (!(aaType == "HOH")) { + lig->addAtom(*at); + lig->setType(aaType); + } + } else { + lig->addAtom(*at); + lig->setType(aaType); + } + }// AminoAcid + else if (tag == "ATOM ") { + // skip N-terminal ACE groups + if (aaType != "ACE") { + // DEBUG: it would be nice to load also alternative atoms + // skip alternative atoms + if (altAaID != '?' && altAaID != '.') { + if (verbose) + cout << "Warning: Skipping extraneous amino acid entry " + << aaNum << " " << atNum << " " << altAaID << ".\n"; + } else { + aa->setType(aaType); + aa->getSideChain().setType(aaType); + + if (!noHAtoms || isHeavyAtom(at->getCode())) { + if (!inSideChain(*aa, *at)) + aa->addAtom(*at); + else { + aa->getSideChain().addAtom(*at); + } + } + } + } else { + if (verbose) + cout << "Warning: Skipping N-terminal ACE group " + << aaNum << " " << atNum << ".\n"; + } + } + delete at; + return aaNum; +} + +/** + * Core function for CIF file parsing. + * @param prot (Protein&) + */ +void CifLoader::loadProtein(Protein& prot) { + PRINT_NAME; + + vector chainList = getAllChains(); + + if (chainList.size() == 0) { + if (verbose) + cout << "Warning: Missing chain ID in the CIF," + "assuming the same chain for the entire file.\n"; + chainList.push_back(char(' ')); + } + + unsigned int readingModel = model; + bool loadChain = false; + + helixCode = ""; + sheetCode = ""; + + string path = "data/AminoAcidHydrogenData.txt"; + const char* inputFile = getenv("VICTOR_ROOT"); + if (inputFile == NULL) + ERROR("Environment variable VICTOR_ROOT was not found.", exception); + + AminoAcidHydrogen::loadParam(((string) inputFile + path).c_str()); + + for (unsigned int i = 0; i < chainList.size(); i++) { + loadChain = false; + // Load all chains + if (allChains) { + loadChain = true; + } else { + // Load only first chain + if (chain == ' ') { + loadChain = true; + chain = '#'; + } + // Load only selected chain + else if (chainList[i] == chain) { + loadChain = true; + chain = '#'; + } + } + + if (loadChain) { + if (verbose) { + cout << "\nLoading chain: ->" << chainList[i] << "<-\n"; + } + setChain(chainList[i]); + + input.clear(); // reset file error flags + input.seekg(0, ios::beg); + + Spacer* sp = new Spacer(); + LigandSet* ls = new LigandSet(); + + string atomLine; + atomLine = readLine(input); + //output << "atomLine: " << atomLine << endl; + + int aaNum = -100000; // infinite negative + int oldAaNum = -100000; + //int lastAa = -10000; + + AminoAcid* aa = new AminoAcid(); + Ligand* lig = new Ligand(); + + int start, end; + + string name = ""; + string tag = ""; + + // read all lines + do { + // read header entry + if (atomLine.find(cif->getTag("header")) != string::npos + && (name == "")) { + sp->setType(cif->getInlineField(atomLine)); + } + // read helix entry + else if (atomLine.find(cif->getTag("helix")) != string::npos) { + cif->parseGroup("helix", atomLine); + + while (atomLine != "# " && input) { + start = stoiDEF(cif->getGroupField("helix", atomLine, + cif->getGroupColumnNumber("helix", "helix start"))); + end = stoiDEF(cif->getGroupField("helix", atomLine, + cif->getGroupColumnNumber("helix", "helix end"))); + helixData.push_back(pair(start, end)); + + helixCode += cif->getGroupField("helix", atomLine, + cif->getGroupColumnNumber("helix", "helix chain")); + + char s[256]; + input.getline(s, 256); + atomLine.assign(s); + //output << "atomLine: " << atomLine << endl; + } + } + // read sheet entry + else if (atomLine.find(cif->getTag("sheet range")) != string::npos) { + cif->parseGroup("sheet range", atomLine); + + while (atomLine != "# " && input) { + start = stoiDEF(cif->getGroupField("sheet range", atomLine, + cif->getGroupColumnNumber("sheet range", "sheet start"))); + end = stoiDEF(cif->getGroupField("sheet range", atomLine, + cif->getGroupColumnNumber("sheet range", "sheet end"))); + sheetData.push_back(pair(start, end)); + + sheetCode += cif->getGroupField("sheet range", atomLine, + cif->getGroupColumnNumber("sheet range", "sheet chain")); + + char s[256]; + input.getline(s, 256); + atomLine.assign(s); + //output << "atomLine: " << atomLine << endl; + } + } + // Parse one line of the "ATOM" and "HETATM" fields + else if (atomLine.substr(0, 6) == "ATOM " || + atomLine.substr(0, 6) == "HETATM") { + tag = atomLine.substr(0, 6); + + // Control model number + readingModel = stouiDEF(cif->getGroupField("atom", atomLine, + cif->getGroupColumnNumber("atom", "model"))); + if (readingModel > model) + break; + // Get only the first model if not specified + if (model == 999) { + model = readingModel; + } + + char chainID = cif->getGroupField("atom", atomLine, + cif->getGroupColumnNumber("atom", "chain")).c_str()[0]; + + if (chainList[i] == chainID) { + if ((model == 999) || (model == readingModel)) { + aaNum = stoiDEF(cif->getGroupField("atom", atomLine, + cif->getGroupColumnNumber("atom", "residue num"))); + + // Insert the Ligand object into LigandSet + if (aaNum != oldAaNum) { + // Print some indexes for the debug + /* + cout << aa->getType1L() << " offset:" << sp->getStartOffset() << " gaps:" + << sp->sizeGaps() << " sizeAmino:" << sp->sizeAmino() << " maxPdbNum:" + << sp->maxPdbNumber() << " aaNum:" << aaNum + << " oldAaNum:" << oldAaNum << " lastAa:" << lastAa << "\n"; + */ + // Skip the first empty AminoAcid + if ((aa->size() > 0) && (aa->getType1L() != 'X')) { + if (sp->sizeAmino() == 0) { + sp->setStartOffset(oldAaNum - 1); + } else { + // Add gaps + //for (int i = lastAa+1; i < oldAaNum; i++){ + for (int i = sp->maxPdbNumber() + 1; i < oldAaNum; i++) { + sp->addGap(i); + } + } + sp->insertComponent(aa); + } + // Ligand + if (lig->size() > 0) { + if (onlyMetalHetAtoms) { + if (lig->isSimpleMetalIon()) { // skip not metal ions + ls->insertComponent(lig); + } + } else { + ls->insertComponent(lig); + } + } + aa = new AminoAcid(); + lig = new Ligand(); + } + oldAaNum = parseCifline(atomLine, tag, lig, aa); + } // end model check + } // end chain check + } + atomLine = readLine(input); + //output << "atomLine: " << atomLine << endl; + } while (input); + + /* + // Print some indexes for the debug + cout << aa->getType1L() << " offset:" << sp->getStartOffset() << " gaps:" + << sp->sizeGaps() << " sizeAmino:" << sp->sizeAmino() << " maxPdbNum:" + << sp->maxPdbNumber() << " aaNum:" << aaNum + << " oldAaNum:" << oldAaNum << " lastAa:" << lastAa << "\n"; + */ + + // last residue/ligand + // AminoAcid + if ((aa->size() > 0) && (aa->getType1L() != 'X')) { + if (sp->sizeAmino() == 0) { + sp->setStartOffset(oldAaNum - 1); + } else { + // Add gaps + //for (int i = lastAa+1; i < oldAaNum; i++){ + for (int i = sp->maxPdbNumber() + 1; i < oldAaNum; i++) { + sp->addGap(i); + } + } + sp->insertComponent(aa); + } + // Ligand + if (lig->size() > 0) { + if (onlyMetalHetAtoms) { + if (lig->isSimpleMetalIon()) { // skip not metal ions + ls->insertComponent(lig); + } + } else { + ls->insertComponent(lig); + } + } + if (verbose) { + cout << "Parsing done\n"; + } + //////////////////////////////////////////////////////////////////// + // Spacer processing + if (sp->sizeAmino() > 0) { + // correct ''fuzzy'' (i.e. incomplete) residues + for (unsigned int j = 0; j < sp->sizeAmino(); j++) { + if ((!sp->getAmino(j).isMember(O)) || + (!sp->getAmino(j).isMember(C)) || + (!sp->getAmino(j).isMember(CA)) || + (!sp->getAmino(j).isMember(N))) { + + // remove residue + sp->deleteComponent(&(sp->getAmino(j))); + + // Add a gap for removed residues + sp->addGap(sp->getStartOffset() + j + 1); + + if (verbose) { + cout << "Warning: Residue number " + << sp->getPdbNumberFromIndex(j) + << " is incomplete and had to be removed.\n"; + } + } + } + if (verbose) { + cout << "Removed incomplete residues\n"; + } + // connect aminoacids + if (!noConnection) { + if (!setBonds(*sp)) { // connect atoms... + valid = false; + if (verbose) { + cout << "Warning: Fail to connect residues in chain: " + << chainList[i] << ".\n"; + } + } + if (verbose) { + cout << "Connected residues\n"; + } + } + + // correct position of leading N atom + sp->setTrans(sp->getAmino(0)[N].getTrans()); + vgVector3 tmp(0.0, 0.0, 0.0); + sp->getAmino(0)[N].setTrans(tmp); + sp->getAmino(0).adjustLeadingN(); + if (verbose) { + cout << "Fixed leading N atom\n"; + } + // Add H atoms + if (!noHAtoms) { + for (unsigned int j = 0; j < sp->sizeAmino(); j++) { + AminoAcidHydrogen::setHydrogen(&(sp->getAmino(j)), false); // second argument is VERBOSE + } + if (verbose) { + cout << "H assigned\n"; + } + if (!noSecondary) { + sp->setDSSP(false); // argument is VERBOSE + if (verbose) { + cout << "DSSP assigned\n"; + } + } + } + + // assign secondary structure from torsion angles + if (!noSecondary) { + assignSecondary(*sp); + if (verbose) { + cout << "Torsional SS assigned\n"; + } + } + } else { + if (verbose) { + cout << "Warning: No residues in chain: " << chainList[i] << ".\n"; + } + } + + //////////////////////////////////////////////////////////////////// + // Load data into protein object + Polymer* pol = new Polymer(); + pol->insertComponent(sp); + if (verbose) { + cout << "Loaded AminoAcids: " << sp->size() << "\n"; + } + if (!(noHetAtoms)) { + //insertion only if LigandSet is not empty + if (ls->sizeLigand() > 0) { + pol->insertComponent(ls); + if (verbose) { + cout << "Loaded Ligands: " << ls->size() << "\n"; + } + } else { + if (verbose) { + cout << "Warning: No ligands in chain: " << chainList[i] << ".\n"; + } + } + } + + prot.addChain(chainList[i]); + prot.insertComponent(pol); + + } // end loadChain + } // chains iteration + +} \ No newline at end of file diff --git a/Biopool/Sources/CifLoader.h b/Biopool/Sources/CifLoader.h new file mode 100644 index 0000000..ab93bed --- /dev/null +++ b/Biopool/Sources/CifLoader.h @@ -0,0 +1,204 @@ +/* This file is part of Victor. + + Victor is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + Victor 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 General Public License for more details. + + You should have received a copy of the GNU General Public License + along with Victor. If not, see . + */ + +#ifndef CIFLOADER_H +#define CIFLOADER_H + +// Includes: +#include +#include +#include +#include +#include +#include +#include + +#include "CifStructure.h" + + +// Global constants, typedefs, etc. (to avoid): + +namespace Victor { + namespace Biopool { + + /** + * @brief Loads components (Atoms, Groups, Spacer, etc.) in standard CIF format. + */ + class CifLoader : public Loader { + public: + + // CONSTRUCTORS/DESTRUCTOR: + /** + * Constructor. + * @param _input = CIF file object + * @param _output = log file + * @param _permissive = if true, allows loading residues with missing atoms + * @param _noHAtoms = if true, doesn't load Hydrogens + * @param _noHetAtoms = if true, doesn't load het atoms + * @param _noSecondary = if true, doesn't load secondary structure (neither the one calculated from torsional angles nor the DSSP) + * @param _noConnection = if true, doesn't connect residues + * @param _noWater = if true, doesn't load water atoms + * @param _verb = if true, verbose mode + * @param _allChains = if true, loads all chains + * @param _NULL = the name of the chain to be loaded, if not provided only loads the first chain + * @param _onlyMetal = if true, load only metals as ligands + * @param _noNucleotideChains = if true, doesn't load DNA/RNA chains + */ + CifLoader(istream& _input = cin, ostream& output = cout, bool _permissive = false, + bool _noHAtoms = false, bool _noHetAtoms = false, + bool _noSecondary = false, bool _noConnection = false, + bool _noWater = true, bool _verb = true, bool _allChains = false, + string _NULL = "", bool _onlyMetal = false, + bool _noNucleotideChains = true); + + // this class uses the implicit copy operator. + + virtual ~CifLoader(); + + // PREDICATES: + + bool isValid(); + void checkModel(); //to check input values ​​ + void checkAndSetChain(); //chosen by the user + unsigned int getMaxModels(); + vector getAllChains(); + + // MODIFIERS: + + void setPermissive(); + void setNonPermissive(); + void setVerbose(); + void setNoVerbose(); + void setChain(char _ch); + void setModel(unsigned int _mod); + void setAltAtom(char _a); + void setNoHAtoms(); + void setNoHetAtoms(); + void setOnlyMetalHetAtoms(); + void setNoSecondary(); + void setWithSecondary(); + void setNoConnection(); + void setWithConnection(); + void setWater(); + void setAllChains(); + + //virtual void loadSpacer(Spacer& sp); + //virtual void loadLigandSet(LigandSet& l); + + virtual void loadProtein(Protein& prot); + + //virtual void loadNucleotideChainSet(NucleotideChainSet& ns); //new class, new code by Damiano + + protected: + // HELPERS: + bool setBonds(Spacer& sp); + bool inSideChain(const AminoAcid& aa, const Atom& at); + void assignSecondary(Spacer& sp); + int parseCifline(string atomLine, string tag, Ligand* lig, AminoAcid* aa); + + // ATTRIBUTES + private: + istream& input; // input stream + ostream& output; + bool permissive; // + bool valid; // + bool noHAtoms; // + bool noHetAtoms; // hetatms contain water, simpleMetalIons and cofactors + bool onlyMetalHetAtoms; // with this flag we select only 2nd cathegory + bool noSecondary; + bool noConnection; // skip connecting aminoacids + bool noWater; // + bool verbose; + bool allChains; // + char chain; // chain ID to be loaded + unsigned int model; // model number to be loaded + char altAtom; // ID of alternate atoms to be loaded + bool noNucleotideChains; // does not load nucleotide atoms + + string helixCode; // parallel vector of helix data, chain name for each helixData element + string sheetCode; + + vector > helixData; // begin and end of the helix + vector > sheetData; + + CifStructure* cif; + }; + + inline bool CifLoader::isValid() { + return valid; + } + + inline void CifLoader::setPermissive() { + permissive = true; + } + + inline void CifLoader::setNonPermissive() { + permissive = false; + } + + inline void CifLoader::setVerbose() { + verbose = true; + } + + inline void CifLoader::setNoVerbose() { + verbose = false; + } + + inline void CifLoader::setChain(char _ch) { + chain = _ch; + } + + inline void CifLoader::setModel(unsigned int _mod) { + model = _mod; + } + + inline void CifLoader::setAltAtom(char _a) { + altAtom = _a; + } + + inline void CifLoader::setNoHAtoms() { + noHAtoms = true; + } + + inline void CifLoader::setNoHetAtoms() { + noHetAtoms = true; + } + + inline void CifLoader::setNoSecondary() { + noSecondary = true; + } + + inline void CifLoader::setWithSecondary() { + noSecondary = false; + } + + inline void CifLoader::setNoConnection() { + noConnection = true; + } + + inline void CifLoader::setWithConnection() { + noConnection = false; + } + + inline void CifLoader::setAllChains() { + allChains = true; + } + + } // namespace Biopool +} // namespace Victor + +#endif /* CIFLOADER_H */ + diff --git a/Biopool/Sources/CifSaver.cc b/Biopool/Sources/CifSaver.cc new file mode 100644 index 0000000..d58b0da --- /dev/null +++ b/Biopool/Sources/CifSaver.cc @@ -0,0 +1,412 @@ +/* This file is part of Victor. + + Victor is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + Victor 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 General Public License for more details. + + You should have received a copy of the GNU General Public License + along with Victor. If not, see . + */ + +// Includes: +#include +#include + +#include "CifSaver.h" + +// Global constants, typedefs, etc. (to avoid): +using namespace Victor; +using namespace Victor::Biopool; +using namespace std; + +// CONSTRUCTORS/DESTRUCTOR: + +CifSaver::CifSaver(ostream& _output) : +output(_output), writeSeq(true), writeSecStr(true), writeTer(true), +atomOffset(0), aminoOffset(0), ligandOffset(0), chain(' '), +atomGroupPrinted(false) { + cif = new CifStructure(_output); +} + +CifSaver::~CifSaver() { + delete cif; + PRINT_NAME; +} + +// PREDICATES: + +// MODIFIERS: + +/** + * Saves a group in CIF format. + * @param group reference + * @return void + */ +void CifSaver::saveGroup(Group& gr) { + gr.sync(); + + if (!atomGroupPrinted) { + cif->printGroup("atom"); + atomGroupPrinted = true; + } + + for (unsigned int i = 0; i < gr.size(); i++) { + string atName = gr[i].getType(); + + // cosmetics: OXT has to be output after + // the sidechain and therefore goes in saveSpacer + if (atName == "OXT") { + continue; + } + + // Added variable for correcting atom type H (last column in PDBs) + char atomOneLetter; + if (!isdigit(atName[0])) { + atomOneLetter = atName[0]; + } else { + atomOneLetter = atName[1]; + } + + // Added control for size by Damiano Piovesan + // example HG12 + if (!isdigit(atName[0]) && (atName.size() < 4)) { + atName = ' ' + atName; + } + while (atName.size() < 4) { + atName += ' '; + } + + // if fields have default values (0 or X), assigns the CIF unknown value (?) + // or a possibly correct value + char asymId = gr[i].getAsymId(); + string entityId = gr[i].getEntityId(); + int model = gr[i].getModel(); + + if (asymId == 'X') { + asymId = '?'; + } + + if (entityId == "0") { + entityId = "?"; + } + + if (model == 0) { + model = 1; + } + + output << setw(7) << left << "ATOM" << + setw(6) << gr[i].getNumber() << + setw(2) << atomOneLetter << + setw(5) << left << atName << + setw(2) << "." << + setw(4) << gr.getType() << + setw(2) << asymId << + setw(2) << entityId << + setw(4) << aminoOffset << + setw(2) << "?" << + setw(8) << setprecision(3) << gr[i].getCoords().x << + setw(8) << setprecision(3) << gr[i].getCoords().y << + setw(8) << setprecision(3) << gr[i].getCoords().z << + setw(6) << setprecision(2) << gr[i].getOccupancy() << + setw(7) << left << setprecision(2) << gr[i].getBFac() << + setw(2) << "?" << + setw(2) << "?" << + setw(2) << "?" << + setw(2) << "?" << + setw(2) << "?" << + setw(2) << "?" << + setw(4) << aminoOffset << + setw(4) << gr.getType() << + setw(2) << chain << + setw(5) << left << atName << + setw(2) << model << + endl; + + atomOffset = gr[i].getNumber() + 1; + } + + //aminoOffset++; +} + +/** + * Saves a sidechain in CIF format. + * @param sideChain side chain to save + */ +void CifSaver::saveSideChain(SideChain& sc) { + saveGroup(sc); +} + +/** + * Saves an aminoacid in CIF format. + * @param AminoAcid aminoacid to save + */ +void CifSaver::saveAminoAcid(AminoAcid& aa) { + saveGroup(aa); +} + +/** + * Saves a spacer in CIF format. + * @param Spacer spacer to save + */ +void CifSaver::saveSpacer(Spacer& sp) { + PRINT_NAME; + + if (sp.size() > 0) { + unsigned int oldPrec = output.precision(); + ios::fmtflags oldFlags = output.flags(); + output.setf(ios::fixed, ios::floatfield); + + //method of class Component. It checks how deep is the spacer + if (sp.getDepth() == 0) { + if (writeTer) { + output << "data_" << sp.getType() << endl; + output << "# " << endl; + output << cif->getTag("header") << " " << sp.getType() + << " " << endl; + output << "# " << endl; + } + aminoOffset = 0; + atomOffset = sp.getAtomStartOffset(); + } + + if (writeSeq) + writeSeqRes(sp); + if (writeSecStr) + writeSecondary(sp); + + aminoOffset = sp.getStartOffset(); + atomOffset = sp.getAtomStartOffset(); + + //saving is one ammino at a time + for (unsigned int i = 0; i < sp.sizeAmino(); i++) { + aminoOffset++; + while ((sp.isGap(aminoOffset)) && (aminoOffset < sp.maxPdbNumber())) { + aminoOffset++; + } + //cout << i << " " << aminoOffset << "\n"; + sp.getAmino(i).save(*this); + } + + // cosmetics: write OXT after last side chain + if (sp.getAmino(sp.sizeAmino() - 1).isMember(OXT)) { + unsigned int index = sp.sizeAmino() - 1; + + // if fields have default values (0 or X), assigns the CIF unknown value (?) + // or a possibly correct value + char asymId = sp.getAmino(index)[OXT].getAsymId(); + string entityId = sp.getAmino(index)[OXT].getEntityId(); + int model = sp.getAmino(index)[OXT].getModel(); + + if (asymId == 'X') { + asymId = '?'; + } + + if (entityId == "0") { + entityId = "?"; + } + + if (model == 0) { + model = 1; + } + + output << setw(7) << left << "ATOM" << + setw(6) << sp.getAmino(index)[OXT].getNumber() << + setw(2) << "O" << + setw(5) << left << "OXT" << + setw(2) << "." << + setw(4) << sp.getAmino(index).getType() << + setw(2) << asymId << + setw(2) << entityId << + setw(4) << aminoOffset << + setw(2) << "?" << + setw(8) << setprecision(3) << sp.getAmino(index)[OXT].getCoords().x << + setw(8) << setprecision(3) << sp.getAmino(index)[OXT].getCoords().y << + setw(8) << setprecision(3) << sp.getAmino(index)[OXT].getCoords().z << + setw(6) << setprecision(2) << sp.getAmino(index)[OXT].getOccupancy() << + setw(7) << left << setprecision(2) << sp.getAmino(index)[OXT].getBFac() << + setw(2) << "?" << + setw(2) << "?" << + setw(2) << "?" << + setw(2) << "?" << + setw(2) << "?" << + setw(2) << "?" << + setw(4) << aminoOffset << + setw(4) << sp.getAmino(index).getType() << + setw(2) << chain << + setw(5) << "OXT" << + setw(2) << model << + endl; + } + + output.precision(oldPrec); + output.flags(oldFlags); + aminoOffset = 0; //necessary if the's more than one spacer + } + +} + +/** + * Saves a Ligand in CIF format. + * @param Ligand ligand to save + */ +void CifSaver::saveLigand(Ligand& gr) { + gr.sync(); + unsigned int oldPrec = output.precision(); + ios::fmtflags oldFlags = output.flags(); + output.setf(ios::fixed, ios::floatfield); + + string aaType = gr.getType(); + + string tag = "HETATM"; + if (isKnownNucleotide(nucleotideThreeLetterTranslator(aaType))) { + tag = "ATOM "; + } + + //print all HETATM of a ligand + for (unsigned int i = 0; i < gr.size(); i++) { + string atType = gr[i].getType(); + aaType = gr.getType(); + string atTypeShort; //last column in a Pdb File + + if (atType != aaType) { + atTypeShort = atType[0]; + } else { + atTypeShort = atType; + } + + // if fields have default values (0 or X), assigns the CIF unknown value (?) + // or a possibly correct value + char asymId = gr[i].getAsymId(); + string entityId = gr[i].getEntityId(); + int model = gr[i].getModel(); + + if (asymId == 'X') { + asymId = '?'; + } + + if (entityId == "0") { + entityId = "?"; + } + + if (model == 0) { + model = 1; + } + + output << setw(7) << left << tag << + setw(6) << gr[i].getNumber() << + setw(2) << atTypeShort << + setw(5) << left << atType << + setw(2) << "." << + setw(4) << aaType << + setw(2) << asymId << + setw(2) << entityId << + setw(4) << ligandOffset << + setw(2) << "?" << + setw(8) << setprecision(3) << gr[i].getCoords().x << + setw(8) << setprecision(3) << gr[i].getCoords().y << + setw(8) << setprecision(3) << gr[i].getCoords().z << + setw(6) << setprecision(2) << gr[i].getOccupancy() << + setw(7) << left << setprecision(2) << gr[i].getBFac() << + setw(2) << "?" << + setw(2) << "?" << + setw(2) << "?" << + setw(2) << "?" << + setw(2) << "?" << + setw(2) << "?" << + setw(4) << aminoOffset << + setw(4) << aaType << + setw(2) << chain << + setw(5) << atType << + setw(2) << model << + endl; + } + output << "# " << endl; + + ligandOffset++; + output.precision(oldPrec); + output.flags(oldFlags); +} + +/** + * Saves a LigandSet in CIF format. + * @param LigandSet set of ligands to save + */ +void CifSaver::saveLigandSet(LigandSet& ls) { + ligandOffset = ls.getStartOffset(); //set the offset for current LigandSet + + for (unsigned int i = 0; i < ls.sizeLigand(); i++) { + while ((ls.isGap(ligandOffset)) + && (ligandOffset < ls.maxPdbNumber())) + ligandOffset++; + ls[i].save(*this); + } +} + +/** + * Saves a Protein in PDB format. + * @param Protein protein to save + */ +void CifSaver::saveProtein(Protein& prot) { + //if (prot.sizeProtein()==0) + // ERROR("Empty Protein",exception); + + Spacer* sp = NULL; + LigandSet* ls = NULL; + + for (unsigned int i = 0; i < prot.sizeProtein(); i++) { + setChain(prot.getChainLetter(i)); //set the actual chain's ID + sp = prot.getSpacer(i); + saveSpacer(*sp); + } + + for (unsigned int i = 0; i < prot.sizeProtein(); i++) { + setChain(prot.getChainLetter(i)); //set the actual chain's ID + ls = prot.getLigandSet(i); + + if (ls != NULL) { + saveLigandSet(*ls); + } + } +} + +/** + * Writes the SEQRES entry (CIF format) for a spacer. + * @param Spacer spacer to write + */ +void CifSaver::writeSeqRes(Spacer& sp) { + cif->printGroup("entity poly"); + + + for (unsigned int i = 0; i < sp.sizeAmino(); i++) { + // if fields have default values (0 or X), assigns the CIF unknown value (?) + // or a possibly correct value + string entityId = sp.getAmino(i).getAtom(0).getEntityId(); + + if (entityId == "0") { + entityId = "?"; + } + + output << setw(2) << left << entityId << + setw(4) << i + 1 << + setw(4) << sp.getAmino(i).getType() << + setw(2) << "n" << + endl; + } + output << "# " << endl; +} + +/** + * Writes the secondary information (PDB format) for a spacer, e.g. HELIX, + * SHEET, etc. + * @param sideChain reference + * @return void + */ +void CifSaver::writeSecondary(Spacer& sp) { + +} \ No newline at end of file diff --git a/Biopool/Sources/CifSaver.h b/Biopool/Sources/CifSaver.h new file mode 100644 index 0000000..49e28cf --- /dev/null +++ b/Biopool/Sources/CifSaver.h @@ -0,0 +1,152 @@ +/* This file is part of Victor. + + Victor is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + Victor 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 General Public License for more details. + + You should have received a copy of the GNU General Public License + along with Victor. If not, see . + */ + +#ifndef CIFSAVER_H +#define CIFSAVER_H + +// Includes: +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "CifStructure.h" + +// Global constants, typedefs, etc. (to avoid): + +namespace Victor { + namespace Biopool { + + /** + * @brief Saves components (Atoms, Groups, etc.) in standard PDB format. + * */ + class CifSaver : public Saver { + public: + // CONSTRUCTORS/DESTRUCTOR: + + /** + * Basic constructor. By default it writes sequence, + * secondary structure and the term line. + * @param _output the output file object + */ + CifSaver(ostream& _output = cout); + + // this class uses the implicit copy operator. + + virtual ~CifSaver(); + + // PREDICATES: + + void endFile(); + + // MODIFIERS: + + void setWriteSecondaryStructure(); + void setDoNotWriteSecondaryStructure(); + void setWriteSeqRes(); + void setDoNotWriteSeqRes(); + void setWriteAtomOnly(); + void setWriteAll(); + void setChain(char _ch); + + /** + * Saves a group in CIF format. + * @param group reference + * @return void + */ + virtual void saveGroup(Group& gr); + virtual void saveSideChain(SideChain& sc); + virtual void saveAminoAcid(AminoAcid& aa); + virtual void saveSpacer(Spacer& sp); + virtual void saveLigand(Ligand& l); + virtual void saveLigandSet(LigandSet& l); + virtual void saveProtein(Protein& prot); + + protected: + + private: + // HELPERS: + + // writes SEQRES entry + void writeSeqRes(Spacer& sp); + + // writes secondary entries (SHEET, HELIX, etc.) + void writeSecondary(Spacer& sp); + + // ATTRIBUTES + ostream& output; // output stream + bool writeSeq, writeSecStr, writeTer; + + // offsets that determine at which atom, + // aminoacid and ligand number to start + unsigned int atomOffset, ligandOffset; + int aminoOffset; + char chain; // chain ID + + bool atomGroupPrinted; + + CifStructure* cif; + }; + + inline void CifSaver::endFile() { + output << "END\n"; + } + + inline void CifSaver::setWriteSecondaryStructure() { + writeSecStr = true; + } + + inline void CifSaver::setDoNotWriteSecondaryStructure() { + writeSecStr = false; + } + + inline void CifSaver::setWriteSeqRes() { + writeSeq = true; + } + + inline void CifSaver::setDoNotWriteSeqRes() { + writeSeq = false; + } + + inline void CifSaver::setWriteAtomOnly() { + writeSecStr = false; + writeSeq = false; + writeTer = false; + } + + inline void CifSaver::setWriteAll() { + writeSecStr = true; + writeSeq = true; + writeTer = true; + } + + inline void CifSaver::setChain(char _ch) { + chain = _ch; + } + + } +} //namespace + +#endif /* CIFSAVER_H */ + diff --git a/Biopool/Sources/CifStructure.cc b/Biopool/Sources/CifStructure.cc new file mode 100644 index 0000000..81daebb --- /dev/null +++ b/Biopool/Sources/CifStructure.cc @@ -0,0 +1,385 @@ +/* This file is part of Victor. + + Victor is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + Victor 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 General Public License for more details. + + You should have received a copy of the GNU General Public License + along with Victor. If not, see . + */ + +#include +#include + +#include +#include + +#include "CifStructure.h" + +using namespace Victor; +using namespace Victor::Biopool; +using namespace std; + +CifStructure::CifStructure(istream& input, ostream& output) : +input(input), output(output) { + initData(); +} + +CifStructure::CifStructure(ostream& output) : +output(output), input(cin) { + initData(); +} + +CifStructure::~CifStructure() { +} + +void CifStructure::initData() { + header = "_entry.id"; + + // atom group + atom = "_atom_site."; + atomId = "id "; + chain = "auth_asym_id "; + asymId = "label_asym_id "; + entityId = "label_entity_id "; + residueIns = "pdbx_PDB_ins_code "; + x = "Cartn_x "; + y = "Cartn_y "; + z = "Cartn_z "; + occupancy = "occupancy "; + tempFactor = "B_iso_or_equiv "; + residueNum = "auth_seq_id "; + residueName = "auth_comp_id "; + atomName = "auth_atom_id "; + model = "pdbx_PDB_model_num "; + + // helix group + helix = "_struct_conf."; + helixStart = "beg_auth_seq_id "; + helixEnd = "end_auth_seq_id "; + helixChainId = "beg_auth_asym_id "; + + // sheet group + sheet = "_struct_sheet."; + sheetOrder = "_struct_sheet_order."; + sheetRange = "_struct_sheet_range."; + sheetHbond = "_pdbx_struct_sheet_hbond."; + sheetStart = "beg_auth_seq_id "; + sheetEnd = "end_auth_seq_id "; + sheetChainId = "beg_auth_asym_id "; + + // flags + atomGroupParsed = false; + helixGroupParsed = false; + sheetGroupParsed = false; + sheetOrderGroupParsed = false; + sheetRangeGroupParsed = false; + sheetHboundgroupParsed = false; +} + +vector& CifStructure::getGroup(string name) { + if (name == "atom") { + return atomGroup; + } else if (name == "helix") { + return helixGroup; + } else if (name == "sheet") { + return sheetGroup; + } else if (name == "sheet order") { + return sheetOrderGroup; + } else if (name == "sheet range") { + return sheetRangeGroup; + } else if (name == "sheet hbond") { + return sheetHbondGroup; + } else { + ERROR("getGroup (CifStructure): invalid group name", exception); + } +} + +string CifStructure::getTag(string name) { + if (name == "header") { + return header; + }// atom group + else if (name == "atom") { + return atom; + } else if (name == "atom id") { + return atomId; + } else if (name == "chain") { + return chain; + } else if (name == "atom asym") { + return asymId; + } else if (name == "atom entity") { + return entityId; + } else if (name == "residue ins") { + return residueIns; + } else if (name == "x") { + return x; + } else if (name == "y") { + return y; + } else if (name == "z") { + return z; + } else if (name == "occupancy") { + return occupancy; + } else if (name == "bfac") { + return tempFactor; + } else if (name == "residue num") { + return residueNum; + } else if (name == "residue name") { + return residueName; + } else if (name == "atom name") { + return atomName; + } else if (name == "model") { + return model; + }// helix group + else if (name == "helix") { + return helix; + } else if (name == "helix start") { + return helixStart; + } else if (name == "helix end") { + return helixEnd; + } else if (name == "helix chain") { + return helixChainId; + }// sheet group + else if (name == "sheet") { + return sheet; + } else if (name == "sheet order") { + return sheetOrder; + } else if (name == "sheet range") { + return sheetRange; + } else if (name == "sheet hbond") { + return sheetHbond; + } else if (name == "sheet start") { + return sheetStart; + } else if (name == "sheet end") { + return sheetEnd; + } else if (name == "sheet chain") { + return sheetChainId; + } else { + ERROR("getTag (CifStructure): invalid tag name", exception); + } +} + +int CifStructure::getGroupColumnNumber(string name, string field) { + int col = -1; + vector& group = getGroup(name); + string tag = getTag(field); + for (vector::iterator it = group.begin(); it != group.end(); it++) { + if (*it == tag) { + col = it - group.begin(); + break; + } + } + return col; +} + +string CifStructure::getGroupField(string name, string& line, int columnNum) { + istringstream iss(line); + vector& group = getGroup(name); + vector fields; + string field; + for (unsigned int i = 0; i < group.size(); i++) { + iss >> field; + fields.push_back(field); + } + if (columnNum != -1) { + return fields[columnNum]; + } else { + return "?"; + } +} + +string CifStructure::getInlineField(string& line) { + istringstream iss(line); + string tag, field; + iss >> tag >> field; + return field; +} + +void CifStructure::parseGroup(string name, string& line) { + bool found = false; + vector& group = getGroup(name); + + // exit the function if the group name is already parsed + if (!isGroupParsed(name)) { + while (input) { + string groupName(getTag(name)); + size_t pos = line.find(groupName); + if (pos != string::npos) { + group.push_back(line.substr(pos + groupName.size(), + line.size() - groupName.size())); + found = true; + } else { + found = false; + } + + // exit the loop when the research of the fields is completed + if (!found && group.size() > 1) { + setParsedFlag(name); + break; + } + + line = readLine(input); + } + } +} + +void CifStructure::setParsedFlag(string name) { + if (name == "atom") { + atomGroupParsed = true; + } else if (name == "helix") { + helixGroupParsed = true; + } else if (name == "sheet") { + sheetGroupParsed = true; + } else if (name == "sheet hbond") { + sheetHboundgroupParsed = true; + } else if (name == "sheet order") { + sheetOrderGroupParsed = true; + } else if (name == "sheet range") { + sheetRangeGroupParsed = true; + } else { + ERROR("setParsedFlag (CifStructure): invalid flag name", exception); + } +} + +bool CifStructure::isGroupParsed(string name) { + if (name == "atom") { + return atomGroupParsed; + } else if (name == "helix") { + return helixGroupParsed; + } else if (name == "sheet") { + return sheetGroupParsed; + } else if (name == "sheet hbond") { + return sheetHboundgroupParsed; + } else if (name == "sheet order") { + return sheetOrderGroupParsed; + } else if (name == "sheet range") { + return sheetRangeGroupParsed; + } else { + ERROR("isGroupParsed (CifStructure): invalid group name", exception); + } +} + +void CifStructure::printGroup(string name) { + + output << "loop_" << endl; + + if (name == "atom") { + output << "_atom_site.group_PDB " << endl << + "_atom_site.id " << endl << + "_atom_site.type_symbol " << endl << + "_atom_site.label_atom_id " << endl << + "_atom_site.label_alt_id " << endl << + "_atom_site.label_comp_id " << endl << + "_atom_site.label_asym_id " << endl << + "_atom_site.label_entity_id " << endl << + "_atom_site.label_seq_id " << endl << + "_atom_site.pdbx_PDB_ins_code " << endl << + "_atom_site.Cartn_x " << endl << + "_atom_site.Cartn_y " << endl << + "_atom_site.Cartn_z " << endl << + "_atom_site.occupancy " << endl << + "_atom_site.B_iso_or_equiv " << endl << + "_atom_site.Cartn_x_esd " << endl << + "_atom_site.Cartn_y_esd " << endl << + "_atom_site.Cartn_z_esd " << endl << + "_atom_site.occupancy_esd " << endl << + "_atom_site.B_iso_or_equiv_esd " << endl << + "_atom_site.pdbx_formal_charge " << endl << + "_atom_site.auth_seq_id " << endl << + "_atom_site.auth_comp_id " << endl << + "_atom_site.auth_asym_id " << endl << + "_atom_site.auth_atom_id " << endl << + "_atom_site.pdbx_PDB_model_num " << endl; + } else if (name == "helix") { + output << "_struct_conf.conf_type_id " << endl << + "_struct_conf.id " << endl << + "_struct_conf.pdbx_PDB_helix_id " << endl << + "_struct_conf.beg_label_comp_id " << endl << + "_struct_conf.beg_label_asym_id " << endl << + "_struct_conf.beg_label_seq_id " << endl << + "_struct_conf.pdbx_beg_PDB_ins_code " << endl << + "_struct_conf.end_label_comp_id " << endl << + "_struct_conf.end_label_asym_id " << endl << + "_struct_conf.end_label_seq_id " << endl << + "_struct_conf.pdbx_end_PDB_ins_code " << endl << + "_struct_conf.beg_auth_comp_id " << endl << + "_struct_conf.beg_auth_asym_id " << endl << + "_struct_conf.beg_auth_seq_id " << endl << + "_struct_conf.end_auth_comp_id " << endl << + "_struct_conf.end_auth_asym_id " << endl << + "_struct_conf.end_auth_seq_id " << endl << + "_struct_conf.pdbx_PDB_helix_class " << endl << + "_struct_conf.details " << endl << + "_struct_conf.pdbx_PDB_helix_length " << endl; + } else if (name == "sheet") { + output << "_struct_sheet.id " << endl << + "_struct_sheet.type " << endl << + "_struct_sheet.number_strands " << endl << + "_struct_sheet.details " << endl; + } else if (name == "sheet range") { + output << "_struct_sheet_range.sheet_id " << endl << + "_struct_sheet_range.id " << endl << + "_struct_sheet_range.beg_label_comp_id " << endl << + "_struct_sheet_range.beg_label_asym_id " << endl << + "_struct_sheet_range.beg_label_seq_id " << endl << + "_struct_sheet_range.pdbx_beg_PDB_ins_code " << endl << + "_struct_sheet_range.end_label_comp_id " << endl << + "_struct_sheet_range.end_label_asym_id " << endl << + "_struct_sheet_range.end_label_seq_id " << endl << + "_struct_sheet_range.pdbx_end_PDB_ins_code " << endl << + "_struct_sheet_range.symmetry " << endl << + "_struct_sheet_range.beg_auth_comp_id " << endl << + "_struct_sheet_range.beg_auth_asym_id " << endl << + "_struct_sheet_range.beg_auth_seq_id " << endl << + "_struct_sheet_range.end_auth_comp_id " << endl << + "_struct_sheet_range.end_auth_asym_id " << endl << + "_struct_sheet_range.end_auth_seq_id " << endl; + } else if (name == "sheet hbond") { + output << "_pdbx_struct_sheet_hbond.sheet_id " << endl << + "_pdbx_struct_sheet_hbond.range_id_1 " << endl << + "_pdbx_struct_sheet_hbond.range_id_2 " << endl << + "_pdbx_struct_sheet_hbond.range_1_label_atom_id " << endl << + "_pdbx_struct_sheet_hbond.range_1_label_comp_id " << endl << + "_pdbx_struct_sheet_hbond.range_1_label_asym_id " << endl << + "_pdbx_struct_sheet_hbond.range_1_label_seq_id " << endl << + "_pdbx_struct_sheet_hbond.range_1_PDB_ins_code " << endl << + "_pdbx_struct_sheet_hbond.range_1_auth_atom_id " << endl << + "_pdbx_struct_sheet_hbond.range_1_auth_comp_id " << endl << + "_pdbx_struct_sheet_hbond.range_1_auth_asym_id " << endl << + "_pdbx_struct_sheet_hbond.range_1_auth_seq_id " << endl << + "_pdbx_struct_sheet_hbond.range_2_label_atom_id " << endl << + "_pdbx_struct_sheet_hbond.range_2_label_comp_id " << endl << + "_pdbx_struct_sheet_hbond.range_2_label_asym_id " << endl << + "_pdbx_struct_sheet_hbond.range_2_label_seq_id " << endl << + "_pdbx_struct_sheet_hbond.range_2_PDB_ins_code " << endl << + "_pdbx_struct_sheet_hbond.range_2_auth_atom_id " << endl << + "_pdbx_struct_sheet_hbond.range_2_auth_comp_id " << endl << + "_pdbx_struct_sheet_hbond.range_2_auth_asym_id " << endl << + "_pdbx_struct_sheet_hbond.range_2_auth_seq_id " << endl; + } else if (name == "entity poly") { + output << "_entity_poly_seq.entity_id " << endl << + "_entity_poly_seq.num " << endl << + "_entity_poly_seq.mon_id " << endl << + "_entity_poly_seq.hetero " << endl; + } else if (name == "pdbx poly") { + output << "_pdbx_poly_seq_scheme.asym_id " << endl << + "_pdbx_poly_seq_scheme.entity_id " << endl << + "_pdbx_poly_seq_scheme.seq_id " << endl << + "_pdbx_poly_seq_scheme.mon_id " << endl << + "_pdbx_poly_seq_scheme.ndb_seq_num " << endl << + "_pdbx_poly_seq_scheme.pdb_seq_num " << endl << + "_pdbx_poly_seq_scheme.auth_seq_num " << endl << + "_pdbx_poly_seq_scheme.pdb_mon_id " << endl << + "_pdbx_poly_seq_scheme.auth_mon_id " << endl << + "_pdbx_poly_seq_scheme.pdb_strand_id " << endl << + "_pdbx_poly_seq_scheme.pdb_ins_code " << endl << + "_pdbx_poly_seq_scheme.hetero " << endl; + } else { + ERROR("printGroup (CifStructure): invalid group name", exception); + } +} diff --git a/Biopool/Sources/CifStructure.h b/Biopool/Sources/CifStructure.h new file mode 100644 index 0000000..cd2baea --- /dev/null +++ b/Biopool/Sources/CifStructure.h @@ -0,0 +1,198 @@ +/* This file is part of Victor. + + Victor is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + Victor 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 General Public License for more details. + + You should have received a copy of the GNU General Public License + along with Victor. If not, see . + */ + +#ifndef CIFSTRUCTURE_H +#define CIFSTRUCTURE_H + + +// Includes: +#include +#include +#include + +using std::string; +using std::istream; +using std::ostream; +using std::vector; + +// Global constants, typedefs, etc. (to avoid): + +namespace Victor { + namespace Biopool { + + /** + * @brief Helper class used to hold information from CIF file. + */ + class CifStructure { + public: + /** + * Constructor. + * @param input input file stream + * @param output output file stream + */ + CifStructure(istream& input, ostream& output = cout); + + /** + * Constructor. + * @param output output file stream + * @param input input file stream + */ + CifStructure(ostream& output); + + /** + * Destructor. + */ + virtual ~CifStructure(); + + /** + * Returns the correct collection by group name. + * @param name name of the CIF group + * @return reference to the collection + */ + vector& getGroup(string name); + + /** + * Returns the tag by name. + * @param name name of tag + * @return CIF tag + */ + string getTag(string name); + + /** + * Returns the column number of the field. + * @param name name of the group + * @param field name of the field + * @return field column number + */ + int getGroupColumnNumber(string name, string field); + + /** + * Returns the field of the line at the columnNum column. + * @param name name of the group + * @param line line of the CIF + * @param columnNum number of column + * @return field at columnNum column + */ + string getGroupField(string name, string& line, int columnNum); + + /** + * Returns the field present in the CIF line. + * @param line + * @return + */ + string getInlineField(string& line); + + /** + * Parses group of CIF fields and creates a vector with columns positions. + * @param name name of the group + */ + void parseGroup(string group, string& line); + + /** + * Sets the flag of the parsed group. + * @param name name of the group + */ + void setParsedFlag(string name); + + /** + * Returns true if the group name is parsed, false otherwise. + * @param name name of the group + * @return true if group is parsed, false otherwise + */ + bool isGroupParsed(string name); + + /** + * Prints group records names into output stream. + * @param name name of the group + */ + void printGroup(string name); + + /** + * Returns the input stream. + * @return input stream + */ + istream& getInput(); + + private: + /** + * Initializes data members. + */ + void initData(); + + // CIF file + istream& input; + ostream& output; + + // CIF tags + string header; + + // atom group + string atom; + string atomId; + string chain; + string asymId; + string entityId; + string residueIns; + string x; + string y; + string z; + string occupancy; + string tempFactor; + string residueNum; + string residueName; + string atomName; + string model; + + // helix group + string helix; + string helixStart; + string helixEnd; + string helixChainId; + + // sheet group + string sheet; + string sheetOrder; + string sheetRange; + string sheetHbond; + string sheetStart; + string sheetEnd; + string sheetChainId; + + // collections of CIF group fields + vector atomGroup; + vector helixGroup; + vector sheetGroup; + vector sheetOrderGroup; + vector sheetRangeGroup; + vector sheetHbondGroup; + + // flags + bool atomGroupParsed; + bool helixGroupParsed; + bool sheetGroupParsed; + bool sheetOrderGroupParsed; + bool sheetRangeGroupParsed; + bool sheetHboundgroupParsed; + }; + + inline istream& CifStructure::getInput() { + return input; + } + } // namespace Biopool +} // namespace Victor + +#endif /* CIFSTRUCTURE_H */ + diff --git a/Biopool/Sources/Makefile b/Biopool/Sources/Makefile index 86b43b9..8a96a1c 100644 --- a/Biopool/Sources/Makefile +++ b/Biopool/Sources/Makefile @@ -32,7 +32,8 @@ SOURCES = Identity.cc SimpleBond.cc Bond.cc \ AminoAcid.cc Spacer.cc IntSaver.cc IntLoader.cc SeqSaver.cc PdbLoader.cc \ PdbSaver.cc SeqLoader.cc IntCoordConverter.cc SeqConstructor.cc Ligand.cc \ LigandSet.cc SolvExpos.cc AminoAcidHydrogen.cc Nucleotide.cc \ - RelLoader.cc XyzSaver.cc RelSaver.cc XyzLoader.cc + RelLoader.cc XyzSaver.cc RelSaver.cc XyzLoader.cc \ + CifStructure.cc CifLoader.cc CifSaver.cc OBJECTS = Identity.o SimpleBond.o Bond.o \ @@ -41,7 +42,8 @@ OBJECTS = Identity.o SimpleBond.o Bond.o \ SeqSaver.o PdbLoader.o PdbSaver.o SeqLoader.o \ IntCoordConverter.o SeqConstructor.o Ligand.o LigandSet.o \ SolvExpos.o Protein.o AminoAcidHydrogen.o Nucleotide.o \ - RelLoader.o XyzSaver.o RelSaver.o XyzLoader.o + RelLoader.o XyzSaver.o RelSaver.o XyzLoader.o \ + CifStructure.o CifLoader.o CifSaver.o TARGETS = diff --git a/Biopool/Sources/PdbSaver.cc b/Biopool/Sources/PdbSaver.cc index 861ff01..cd1810d 100644 --- a/Biopool/Sources/PdbSaver.cc +++ b/Biopool/Sources/PdbSaver.cc @@ -274,15 +274,17 @@ void PdbSaver::saveProtein(Protein& prot) { */ void PdbSaver::writeSeqRes(Spacer& sp) { for (unsigned int i = 0; i < sp.sizeAmino() / 13; i++) { - output << "SEQRES " << setw(3) << i << " " << setw(3) - << sp.sizeAmino() << " "; + output << "SEQRES " << setw(3) << i << " " << + setw(1) << chain << " " << + setw(4) << sp.sizeAmino() << " "; for (unsigned int j = 0; j < 13; j++) output << sp.getAmino((i * 13) + j).getType() << " "; output << "\n"; } if (sp.sizeAmino() % 13 > 0) { - output << "SEQRES " << setw(3) << sp.sizeAmino() / 13 + 1 << " " - << setw(3) << sp.sizeAmino() << " "; + output << "SEQRES " << setw(3) << sp.sizeAmino() / 13 + 1 << " " << + setw(1) << chain << " " << + setw(4) << sp.sizeAmino() << " "; for (unsigned int j = 13 * (sp.sizeAmino() / 13); j < sp.sizeAmino(); j++) output << sp.getAmino(j).getType() << " "; output << "\n"; diff --git a/Biopool/Tests/Makefile b/Biopool/Tests/Makefile index 6f35654..d740be5 100644 --- a/Biopool/Tests/Makefile +++ b/Biopool/Tests/Makefile @@ -26,15 +26,20 @@ INC_PATH = -I. # Objects and headers # -SOURCES = TestBiopool.cc TestAtom.h TestAminoAcid.h TestGroup.h TestSpacer.h +#SOURCES = TestBiopool.cc TestAtom.h TestAminoAcid.h TestGroup.h TestSpacer.h +SOURCES = TestCif.cc TestCifLoader.h TestCifStructure.h -OBJECTS = $(SOURCES:.cpp=.o) +#OBJECTS = $(SOURCES:.cpp=.o) +OBJECTS = $(SOURCES:.cc=.o) -TARGETS = TestBiopool +#TARGETS = TestBiopool +TARGETS = TestCif -EXECS = TestBiopool +#EXECS = TestBiopool +EXECS = TestCif -LIBRARY = TESTlibBiopool.a +#LIBRARY = TESTlibBiopool.a +LIBRARY = TESTlibCif.a # # Install rule diff --git a/Biopool/Tests/TestCif.cc b/Biopool/Tests/TestCif.cc new file mode 100644 index 0000000..430ce24 --- /dev/null +++ b/Biopool/Tests/TestCif.cc @@ -0,0 +1,32 @@ +/* + * File: TestCif.cc + * Author: marco + * + * Created on 9-giu-2015, 10.45.40 + */ + +#include +#include + +#include +#include + +using std::cout; +using std::endl; + +int main() { + + CppUnit::TextUi::TestRunner runner; + + cout << "Creating Test Suites:" << endl; + + runner.addTest(TestCifLoader::suite()); + runner.addTest(TestCifStructure::suite()); + + cout << "Running the unit tests." << endl; + + runner.run(); + + return 0; + +} diff --git a/Biopool/Tests/TestCifLoader.h b/Biopool/Tests/TestCifLoader.h new file mode 100644 index 0000000..8ca520c --- /dev/null +++ b/Biopool/Tests/TestCifLoader.h @@ -0,0 +1,96 @@ +/* This file is part of Victor. + + Victor is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + Victor 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 General Public License for more details. + + You should have received a copy of the GNU General Public License + along with Victor. If not, see . + */ + +#ifndef TESTCIFLOADER_H +#define TESTCIFLOADER_H + +#include +#include + +#include +#include +#include +#include +#include + +#include + +using namespace std; +using namespace Victor::Biopool; +using namespace CppUnit; + +class TestCifLoader : public TestFixture { +public: + + TestCifLoader() { + } + + virtual ~TestCifLoader() { + } + + static Test* suite() { + TestSuite* suiteOfTests = new TestSuite("TestCifLoader"); + + suiteOfTests->addTest(new TestCaller("Test 1 - Get max numebers of models", + &TestCifLoader::testGetMaxModels)); + + suiteOfTests->addTest(new TestCaller("Test 2 - Get all chain ids", + &TestCifLoader::testGetAllChains)); + + return suiteOfTests; + } + + void setUp() { + // inizialize CifLoader + string path = getenv("VICTOR_ROOT"); + string input = path + "Biopool/Tests/data/modelTest.cif"; + inFile = new ifstream(input.c_str()); + cl = new CifLoader(*inFile); + + // initialize test parameters + maxModel = 5; + chainIds.push_back('A'); + chainIds.push_back('B'); + chainIds.push_back('C'); + } + + void tearDown() { + delete cl; + delete inFile; + } + +private: + + void testGetMaxModels() { + int max = cl->getMaxModels(); + CPPUNIT_ASSERT_EQUAL(maxModel, max); + } + + void testGetAllChains() { + vector testChain = cl->getAllChains(); + CPPUNIT_ASSERT(chainIds[0] == testChain[0] && + chainIds[1] == testChain[1] && + chainIds[2] == testChain[2]); + } + + int maxModel; + vector chainIds; + + ifstream* inFile; + CifLoader* cl; +}; +#endif /* TESTCIFLOADER_H */ + diff --git a/Biopool/Tests/TestCifStructure.h b/Biopool/Tests/TestCifStructure.h new file mode 100644 index 0000000..507aea2 --- /dev/null +++ b/Biopool/Tests/TestCifStructure.h @@ -0,0 +1,116 @@ +/* This file is part of Victor. + + Victor is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + Victor 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 General Public License for more details. + + You should have received a copy of the GNU General Public License + along with Victor. If not, see . + */ + +#ifndef TESTCIFSTRUCTURE_H +#define TESTCIFSTRUCTURE_H + +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include + +using namespace std; +using namespace Victor::Biopool; +using namespace CppUnit; + +class TestCifStructure : public TestFixture { +public: + + TestCifStructure() { + } + + virtual ~TestCifStructure() { + } + + static Test* suite() { + TestSuite* suiteOfTests = new TestSuite("TestCifLoader"); + + suiteOfTests->addTest(new TestCaller("Get group column number", + &TestCifStructure::testGetGroupColumnNumber)); + + suiteOfTests->addTest(new TestCaller("Get group field", + &TestCifStructure::testGetGroupField)); + + suiteOfTests->addTest(new TestCaller("Get inline field", + &TestCifStructure::testGetInlineField)); + + return suiteOfTests; + } + + void setUp() { + // inizialize CifLoader + string path = getenv("VICTOR_ROOT"); + string input = path + "Biopool/Tests/data/modelTest.cif"; + inFile = new ifstream(input.c_str()); + cs = new CifStructure(*inFile); + + // initialize test parameters + idColumn = 1; + atomId = 1; + headerLine = "_entry.id 25C8 "; + headerId = "25C8"; + } + + void tearDown() { + delete cs; + delete inFile; + } + +private: + + void testGetGroupField() { + string line = readLine(cs->getInput()); + cs->parseGroup("atom", line); + int field = stoiDEF(cs->getGroupField("atom", line, + cs->getGroupColumnNumber("atom", "atom id"))); + + CPPUNIT_ASSERT_EQUAL(atomId, field); + } + + void testGetGroupColumnNumber() { + string line = readLine(cs->getInput()); + cs->parseGroup("atom", line); + int col = cs->getGroupColumnNumber("atom", "atom id"); + + CPPUNIT_ASSERT_EQUAL(idColumn, col); + } + + void testGetInlineField() { + string field = cs->getInlineField(headerLine); + + CPPUNIT_ASSERT_EQUAL(headerId, field); + } + + int idColumn; + int atomId; + string headerLine; + string headerId; + + ifstream* inFile; + CifStructure* cs; + +}; + +#endif /* TESTCIFSTRUCTURE_H */ + diff --git a/Biopool/Tests/data/modelTest.cif b/Biopool/Tests/data/modelTest.cif new file mode 100644 index 0000000..6af1b88 --- /dev/null +++ b/Biopool/Tests/data/modelTest.cif @@ -0,0 +1,43 @@ +loop_ +_atom_site.group_PDB +_atom_site.id +_atom_site.type_symbol +_atom_site.label_atom_id +_atom_site.label_alt_id +_atom_site.label_comp_id +_atom_site.label_asym_id +_atom_site.label_entity_id +_atom_site.label_seq_id +_atom_site.pdbx_PDB_ins_code +_atom_site.Cartn_x +_atom_site.Cartn_y +_atom_site.Cartn_z +_atom_site.occupancy +_atom_site.B_iso_or_equiv +_atom_site.Cartn_x_esd +_atom_site.Cartn_y_esd +_atom_site.Cartn_z_esd +_atom_site.occupancy_esd +_atom_site.B_iso_or_equiv_esd +_atom_site.pdbx_formal_charge +_atom_site.auth_seq_id +_atom_site.auth_comp_id +_atom_site.auth_asym_id +_atom_site.auth_atom_id +_atom_site.pdbx_PDB_model_num +ATOM 1 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? A ? 1 +ATOM 2 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? B ? 1 +ATOM 3 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? C ? 1 +ATOM 4 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? A ? 2 +ATOM 5 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? B ? 2 +ATOM 6 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? C ? 2 +ATOM 7 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? A ? 3 +ATOM 8 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? B ? 3 +ATOM 9 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? C ? 3 +ATOM 10 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? A ? 4 +ATOM 11 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? B ? 4 +ATOM 12 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? C ? 4 +ATOM 13 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? A ? 5 +ATOM 14 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? B ? 5 +ATOM 15 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? C ? 5 +# \ No newline at end of file diff --git a/Makefile.global b/Makefile.global index 1f78a89..0a9ddda 100644 --- a/Makefile.global +++ b/Makefile.global @@ -38,7 +38,7 @@ DEBUGFLAGS = -g NORMALFLAGS = $(DEBUGFLAGS) -FASTFLAGS = -O3 -ffast-math -DNDEBUG -ftemplate-depth-36 +FASTFLAGS = -O3 -ffast-math -DNDEBUG -ftemplate-depth-48 QUICKFLAGS = -DNDEBUG diff --git a/tools/String2Number.cc b/tools/String2Number.cc index ec41c69..68cf368 100644 --- a/tools/String2Number.cc +++ b/tools/String2Number.cc @@ -43,7 +43,7 @@ int stoiDEF(const string &s) { * @Description if the input is no number (e.g. char or string) then error, * if wrong input (integer expected but e.g. float given) then error * changes string into integer: - * */ + */ unsigned int stouiDEF(const string &s) { unsigned int i; double temp; @@ -201,9 +201,8 @@ vector sToVectorOfUIntDEF(const string& s_orig) { } /** - * @Description changes integer into sting: + * @Description changes integer into string: */ - string itosDEF(const int &i) { string s; @@ -213,8 +212,8 @@ string itosDEF(const int &i) { } /** - * @Description changes unsigned integer into sting: - * */ + * @Description changes unsigned integer into string: + */ string uitosDEF(const unsigned int &i) { string s; @@ -224,8 +223,8 @@ string uitosDEF(const unsigned int &i) { } /** - * @Description changes long into sting: - * */ + * @Description changes long into string: + */ string ltosDEF(const long &l) { string s; @@ -235,8 +234,8 @@ string ltosDEF(const long &l) { } /** - * @Description changes float into sting: - * */ + * @Description changes float into string: + */ string ftosDEF(const float &f) { string s; @@ -246,8 +245,8 @@ string ftosDEF(const float &f) { } /** - * @Description changes double into sting: - * */ + * @Description changes double into string: + */ string dtosDEF(const double &d) { string s; @@ -258,7 +257,7 @@ string dtosDEF(const double &d) { /** * @Description tokenize text - * */ + */ vector getTokens(const string& text) { istringstream ist(text.c_str()); char* charLine = new char[text.size() + 1]; // size of string