| #include "ChemMolecule.h"
|
| #include "SASurface.h"
|
|
|
| |
| |
| |
|
|
| const float ChemMolecule::flexProbs[MAX_AA]={0.24, 0.38, 0.14, 0.23, 0.06,
|
| 0.11, 0.05, 0.05, 0.06, 0.07,
|
| 0.05, 0.00, 0.13, 0.13, 0.08,
|
| 0.015,0.02, 0.24, 0.07, 0.02,
|
| 0.0
|
| };
|
|
|
| const float ChemMolecule::maxCaDist[MAX_AA]={8.0, 6.3, 5.5, 5.3, 4.2,
|
| 4.2, 3.0, 2.0, 3.0, 3.0,
|
| 5.4, 3.0, 4.5, 4.4, 3.1,
|
| 5.7, 3.5, 5.7, 7.0, 6.6,
|
| 0.0
|
| };
|
|
|
| bool compareFragment(ChemMolecule::Fragment f1, ChemMolecule::Fragment f2) {
|
| if(f1.first != f2.first) {
|
| return f1.first < f2.first;
|
| }
|
|
|
| return f1.second.first < f2.second.first;
|
| }
|
|
|
| void ChemMolecule::loadMolecule(std::istream &molFile, const ChemLib& chemLib, const PDB::Selector& selector) {
|
| Molecule<ChemAtom>::readAllPDBfile(molFile, selector);
|
| for(Molecule<ChemAtom>::iterator molIter = begin(); molIter != end(); molIter++) {
|
| if(!molIter->isHydrogen()) {
|
| const ChemEntry* entry=chemLib.getLibEntry(molIter->residueName(), molIter->type());
|
| if(entry != NULL) {
|
| molIter->setRadius(entry->getEntryRadius());
|
| molIter->setCharge(entry->getEntryCharge());
|
| molIter->setChemType(entry->getEntryChemType());
|
| }
|
| }
|
| }
|
|
|
| centroid_ = centroid();
|
| setIsProtein();
|
| }
|
|
|
| void ChemMolecule::readAllPDBfile(std::istream &molFile, const PDB::Selector& selector) {
|
| Molecule<ChemAtom>::readAllPDBfile(molFile, selector);
|
| computeASAperAtom();
|
| centroid_ = centroid();
|
| setIsProtein();
|
| }
|
|
|
| float ChemMolecule::computeASAperAtom() {
|
| SASurface sas(*this, 1.8, 10);
|
| float asa = sas.computeASAForAtoms();
|
| Logger::infoMessage() << "ASA of the molecule = " << asa << std::endl;
|
| return asa;
|
| }
|
|
|
| void ChemMolecule::computeSteadyProbs() {
|
| Vector3 caAtom;
|
| for(Molecule<ChemAtom>::iterator molIter = begin(); molIter != end(); molIter++) {
|
| RESIDUE_INDEX resIndex = RESIDUE_INDEX(*molIter);
|
| if(molIter->isBackbone()) {
|
|
|
| if(resIndex == GLY)
|
| molIter->setSteadyProb(1-flexProbs[resIndex]);
|
| else {
|
| molIter->setSteadyProb(1-flexProbs[OTHER]);
|
| if(molIter->isCA())
|
| caAtom = molIter->position();
|
| }
|
| } else {
|
| float distanceFromCa = caAtom.dist(molIter->position());
|
| float normFactor= distanceFromCa/maxCaDist[resIndex];
|
| molIter->setSteadyProb(1-flexProbs[resIndex]*normFactor);
|
| }
|
| }
|
| }
|
|
|
| void ChemMolecule::setIsProtein() {
|
| for(Molecule<ChemAtom>::iterator molIter = begin(); molIter != end(); molIter++) {
|
| if(molIter->isCA()) {
|
| isProtein_ = true;
|
| return;
|
| }
|
| }
|
| isProtein_ = false;
|
| }
|
|
|
| unsigned int ChemMolecule::readBindingSite(const std::string fileName) {
|
| bindingSites_.push_back(std::set<unsigned int>());
|
| return readSiteFile(fileName, bindingSites_.back());
|
| }
|
|
|
| unsigned int ChemMolecule::readBlockingSite(const std::string fileName) {
|
| blockingSites_.push_back(std::set<unsigned int>());
|
| return readSiteFile(fileName, blockingSites_.back());
|
| }
|
|
|
| unsigned int ChemMolecule::readSiteFile(const std::string fileName, std::set<unsigned int>& site) {
|
| std::ifstream file(fileName.c_str());
|
| if(!file) {
|
| std::cerr << "Can't open file " << fileName << std::endl;
|
| Logger::errorMessage() << "Can't open file " << fileName << std::endl;
|
| exit(1);
|
| }
|
|
|
| char chainId;
|
| std::string residueSequenceID;
|
|
|
| while(!file.eof()) {
|
| std::string line;
|
| getline(file, line);
|
| boost::trim(line);
|
| if(line.length() == 0) continue;
|
| int entry = -1;
|
|
|
| char record[line.size()+1];
|
| strcpy(record, line.c_str());
|
|
|
| char* splittedRecord = strtok(record, " ");
|
| residueSequenceID = std::string(splittedRecord);
|
|
|
| splittedRecord = strtok (NULL, " ");
|
| if(splittedRecord != NULL) {
|
| chainId = splittedRecord[0];
|
| } else {
|
| chainId = ' ';
|
| }
|
|
|
| entry = residueEntry(chainId, residueSequenceID);
|
|
|
| if(entry != -1) {
|
| site.insert(entry);
|
| } else {
|
| Logger::errorMessage() << "Invalid line in site file: " << line << std::endl;
|
| std::cerr << "Invalid line in site file: " << line << std::endl;
|
| exit(1);
|
| }
|
| }
|
|
|
| Logger::infoMessage() << site.size() << " residues were read from site file: " << fileName << std::endl;
|
| file.close();
|
| return site.size();
|
| }
|
|
|
|
|
| const std::set<unsigned int>& ChemMolecule::getBindingSite(unsigned int siteNumber) const {
|
| const static std::set<unsigned int> emptySet;
|
| if(bindingSites_.size() > siteNumber) {
|
| return bindingSites_[siteNumber];
|
| }
|
| return emptySet;
|
| }
|
|
|
| const std::set<unsigned int>& ChemMolecule::getBlockingSite(unsigned int siteNumber) const {
|
| const static std::set<unsigned int> emptySet;
|
| if(blockingSites_.size() > siteNumber) {
|
| return blockingSites_[siteNumber];
|
| }
|
| return emptySet;
|
| }
|
|
|
| const ChemAtom& ChemMolecule::getChemAtom(int atomIndex) const {
|
| for(Molecule<ChemAtom>::const_iterator molIter = begin(); molIter != end(); molIter++) {
|
| if((int)molIter->atomIndex() == atomIndex) {
|
| return *molIter;
|
| }
|
| }
|
| Logger::errorMessage() << "Input Error: Can't find atom with atom index " << atomIndex << "!" << std::endl;
|
| std::cerr << "Input Error: Can't find atom with atom index " << atomIndex << "!" << std::endl;
|
| exit(1);
|
| }
|
|
|
| double ChemMolecule::getVolume() const {
|
| double volume = 0.0;
|
| double c = (4.0/3.0)*M_PI;
|
| for(Molecule<ChemAtom>::const_iterator molIter = begin(); molIter != end(); molIter++) {
|
| double r = molIter->getRadius();
|
| volume += c * r * r * r;
|
| }
|
| return volume;
|
| }
|
|
|
| const std::vector<ChemMolecule::Fragment>& ChemMolecule::getFragments() const {
|
| if(fragments_.size() == 0) {
|
| ChemMolecule* nonConstThis = const_cast<ChemMolecule*>(this);
|
| nonConstThis->computeFragments();
|
| }
|
| return fragments_;
|
| }
|
|
|
| void ChemMolecule::computeFragments() {
|
|
|
| char currChain;
|
| int firstResIndex, prevResIndex;
|
| bool currChainSet = false;
|
| for(auto i = begin(); i != end(); i++) {
|
| if(!i->isCA() || i->getAtomEntryType() == HETATM) continue;
|
| char chain = i->chainId();
|
| int resIndex = i->residueIndex();
|
|
|
|
|
| if(currChainSet && currChain == chain &&
|
| resIndex > prevResIndex) {
|
| prevResIndex = resIndex;
|
| } else {
|
| if(currChainSet) {
|
| ResidueRange range(firstResIndex, prevResIndex);
|
| Fragment fragment(std::make_pair(currChain, range));
|
| fragments_.push_back(fragment);
|
| }
|
|
|
| currChain = chain;
|
| firstResIndex = prevResIndex = resIndex;
|
| currChainSet = true;
|
| }
|
| }
|
|
|
| if(currChainSet) {
|
| ResidueRange range(firstResIndex, prevResIndex);
|
| Fragment fragment(std::make_pair(currChain, range));
|
| fragments_.push_back(fragment);
|
| }
|
|
|
| std::sort(fragments_.begin(), fragments_.end(), compareFragment);
|
|
|
| for(int i=0; i<(int)fragments_.size(); i++) {
|
| std::cout << "Fragment " << i << " chainId "
|
| << fragments_[i].first << " range "
|
| << fragments_[i].second.first << ":"
|
| << fragments_[i].second.second << std::endl;
|
| }
|
| }
|
|
|
| void ChemMolecule::addMol2Type() {
|
|
|
|
|
| for(auto i = begin(); i != end(); i++) {
|
| const char* type = i->type();
|
| char res = i->residueType();
|
| i->setMol2Type(UNK_MOL2_TYPE);
|
| if(i->isN()) i->setMol2Type(N2);
|
| if(i->isC()) i->setMol2Type(C2);
|
| if(i->isCA()) i->setMol2Type(C3);
|
| if(i->isO()) i->setMol2Type(O2);
|
| if(i->isCB()) i->setMol2Type(C3);
|
|
|
|
|
| if(i->isCD() && (res == 'K' || res == 'P' || res == 'R')) i->setMol2Type(C3);
|
| if(i->isCD() && (res == 'Q' || res == 'E' )) i->setMol2Type(C2);
|
|
|
| if(i->isCG() && (res == 'K' || res == 'P' || res == 'R' || res == 'M' || res == 'L' || res == 'Q' || res == 'E' )) i->setMol2Type(C3);
|
| if(i->isCG() && (res == 'D' || res == 'N')) i->setMol2Type(C2);
|
| if(i->isCG() && (res == 'H' || res == 'F' || res == 'W' || res == 'Y')) i->setMol2Type(Car);
|
| if(i->isCZ() && (res == 'F' || res == 'Y')) i->setMol2Type(Car);
|
| if(i->isCZ() && res == 'R') i->setMol2Type(Ccat);
|
| if((i->isCD1() || i->isCD2()) && (res == 'F' || res == 'W' || res == 'Y' || res == 'H')) i->setMol2Type(Car);
|
| if((i->isCD1() || i->isCD2()) && (res == 'L' || res == 'I')) i->setMol2Type(C3);
|
| if(i->isCE() || i->isCG1() || i->isCG2()) i->setMol2Type(C3);
|
| if(i->isCE1() || i->isCE2() || i->isCE3() || i->isCZ2() || i->isCZ3() || i->isCH2() ) i->setMol2Type(Car);
|
|
|
|
|
| if (i->isSugarCarbon())
|
| i->setMol2Type(C3);
|
|
|
|
|
| if((i->isOD1() || i->isOE1()) && ((res == 'N' || res == 'Q'))) i->setMol2Type(O2);
|
| if((i->isOD1() || i->isOE1()|| i->isOD2()|| i->isOE2()) && (res == 'D' || res == 'E')) i->setMol2Type(Oco2);
|
| if(i->isOH() || i->isOG()|| i->isOG1()) i->setMol2Type(O3);
|
|
|
|
|
| if (i->isSugarOxygen())
|
| i->setMol2Type(O3);
|
| if (i->isPhosphateOxygen())
|
| i->setMol2Type(Oco2);
|
|
|
|
|
| if(i->isNH1() || i->isNH2() || i->isNE()) i->setMol2Type(Npl3);
|
| if((i->isND2() && res == 'N') || (i->isNE2() && res == 'Q')) i->setMol2Type(Nam);
|
| if(((i->isND1() || i->isNE2()) && res == 'H') || (i->isNE1() && res == 'W')) i->setMol2Type(Nar);
|
| if(i->isNZ()) i->setMol2Type(N4);
|
|
|
|
|
| if (i->isNitrogenousBaseAtom()) {
|
| if (type[1] == 'N') {
|
| if ((type[2] == '4' && res == 'C') || (type[2] == '6' && res == 'A') || (type[2] == '2' && res == 'G'))
|
| i->setMol2Type(Npl3);
|
| else
|
| i->setMol2Type(Nar);
|
| }
|
| else if (type[1] == 'O')
|
| i->setMol2Type(O2);
|
| else if (type[1] == 'C')
|
| i->setMol2Type(Car);
|
| }
|
|
|
|
|
| if(i->isSD() || i->isSG()) i->setMol2Type(S3);
|
|
|
|
|
| if(i->isP()) i->setMol2Type(P3);
|
|
|
| if(i->getMol2Type() == UNK_MOL2_TYPE)
|
| std::cout << "Can't assign mol2 type:" << type[0] << type[1] << type[2] << type[3] << ":" << std::endl;
|
| }
|
| }
|
|
|
