| #include "SuperBB.h"
|
| #include "HierarchicalFold.h"
|
|
|
|
|
| SuperBB::SuperBB(std::shared_ptr<const BB> bb)
|
| : backBonePen_(0), transScore_(0), weightedTransScore_(100) {
|
| bbs_.push_back(bb);
|
| Vector3 v(0, 0, 0);
|
| Matrix3 M(1);
|
| RigidTrans3 T(M, v);
|
| trans_.push_back(T);
|
| size_ = 1;
|
| bitIDS_ = bb->bitId();
|
| }
|
|
|
| float getWeightedTransScore(std::vector<FoldStep> steps, std::vector<std::shared_ptr<const BB>> bbs) {
|
| unsigned int totalCa = 0;
|
| float totalScore = 0;
|
| std::map<unsigned int, unsigned int> bbIdToCaLength;
|
| for (std::shared_ptr<const BB> bb : bbs) {
|
| bbIdToCaLength[bb->getID()] = bb->getNumOfAtoms();
|
| }
|
|
|
| for (FoldStep step : steps) {
|
| std::vector<unsigned int> right;
|
| std::vector<unsigned int> left;
|
|
|
| right.push_back(step.i_);
|
| left.push_back(step.j_);
|
|
|
| std::vector<FoldStep> notUsedSteps(steps);
|
|
|
| while (!notUsedSteps.empty()) {
|
| FoldStep step2 = notUsedSteps.back();
|
| notUsedSteps.pop_back();
|
| if (step2.i_ == step.i_ && step2.j_ == step.j_)
|
| continue;
|
| if (std::find(right.begin(), right.end(), step2.i_) != right.end() &&
|
| std::find(right.begin(), right.end(), step2.j_) == right.end()) {
|
| right.push_back(step2.j_);
|
| } else if (std::find(right.begin(), right.end(), step2.j_) != right.end() &&
|
| std::find(right.begin(), right.end(), step2.i_) == right.end()) {
|
| right.push_back(step2.i_);
|
| } else if (std::find(left.begin(), left.end(), step2.i_) != left.end() &&
|
| std::find(left.begin(), left.end(), step2.j_) == left.end()) {
|
| left.push_back(step2.j_);
|
| } else if (std::find(left.begin(), left.end(), step2.j_) != left.end() &&
|
| std::find(left.begin(), left.end(), step2.i_) == left.end()) {
|
| left.push_back(step2.i_);
|
| } else {
|
| notUsedSteps.insert(notUsedSteps.begin(), step2);
|
| }
|
| }
|
|
|
| unsigned int rightSize = 0;
|
| unsigned int leftSize = 0;
|
|
|
| for (unsigned int elem : right)
|
| rightSize += bbIdToCaLength[elem];
|
|
|
| for (unsigned int elem : left)
|
| leftSize += bbIdToCaLength[elem];
|
|
|
| unsigned int usedSize = rightSize;
|
| if (rightSize > leftSize)
|
| usedSize = leftSize;
|
|
|
| totalScore += usedSize * step.tScore_;
|
| totalCa += usedSize;
|
| }
|
| return totalScore / totalCa;
|
| }
|
|
|
| void SuperBB::join(const RigidTrans3 &trans, const SuperBB &other, int bbPen, FoldStep &step, float transScore) {
|
|
|
| for (unsigned int j = 0; j < other.size_; j++) {
|
| bbs_.push_back(other.bbs_[j]);
|
| trans_.push_back(trans * other.trans_[j]);
|
| }
|
| size_ += other.size_;
|
| bitIDS_ |= other.bitIDS_;
|
| backBonePen_ = bbPen;
|
|
|
| transScore_ += other.transScore_ + transScore;
|
|
|
| if (!other.foldSteps_.empty()) {
|
| foldSteps_.insert(foldSteps_.end(), other.foldSteps_.begin(), other.foldSteps_.end());
|
| }
|
| foldSteps_.push_back(step);
|
|
|
| weightedTransScore_ = getWeightedTransScore(foldSteps_, bbs_);
|
| }
|
|
|
| void SuperBB::replaceIdentBB(BitId oldBBBitId, std::shared_ptr<const BB> bb) {
|
| if((bb->bitId() & bitIDS_) != 0)
|
| throw std::runtime_error("Tried to replace BB with BB that already exists in SuperBB");
|
|
|
| std::shared_ptr<const BB> oldBB;
|
| unsigned int bbIndex = -1;
|
| for (unsigned int i = 0; i < size_; i++) {
|
| if (bbs_[i]->bitId() == oldBBBitId) {
|
| oldBB = bbs_[i];
|
| bbIndex = i;
|
| break;
|
| }
|
| }
|
| if(bbIndex == -1)
|
| throw std::runtime_error("SuperBB::replaceIdentBB: oldBBBitId not found in SuperBB");
|
|
|
| bbs_[bbIndex] = bb;
|
|
|
|
|
| bitIDS_[oldBB->getID()] = false;
|
| bitIDS_[bb->getID()] = true;
|
|
|
| for(FoldStep &step : foldSteps_) {
|
| if(step.i_ == oldBB->getID())
|
| step.i_ = bb->getID();
|
| if(step.j_ == oldBB->getID())
|
| step.j_ = bb->getID();
|
| }
|
| }
|
|
|
| bool SuperBB::isPenetrating(const RigidTrans3 &trans, const SuperBB &other, float threshold) const {
|
| for (unsigned int i = 0; i < size_; i++) {
|
| RigidTrans3 t = (!trans_[i]) * trans;
|
| for (unsigned int j = 0; j < other.size_; j++) {
|
| if (bbs_[i]->isPenetrating(t * other.trans_[j], *other.bbs_[j], threshold))
|
| return true;
|
| }
|
| }
|
| return false;
|
| }
|
|
|
|
|
| double SuperBB::calcRmsd(const SuperBB &other, std::vector<std::vector<unsigned int>> &identGroups) const {
|
| std::vector<std::vector<unsigned int>> presentIdentGroups;
|
|
|
| for (std::vector<unsigned int> identGroup : identGroups) {
|
| std::vector<unsigned int> possiblyRelevantIdentGroup;
|
| for (unsigned int i : identGroup) {
|
| if (bitIDS_.test(i))
|
| possiblyRelevantIdentGroup.push_back(i);
|
| }
|
| if (possiblyRelevantIdentGroup.size() >= 2)
|
| presentIdentGroups.push_back(possiblyRelevantIdentGroup);
|
| }
|
|
|
| std::map<unsigned int, unsigned int> bbIdToThisBBIndex;
|
| std::map<unsigned int, unsigned int> bbIdToOtherBBIndex;
|
| std::vector<unsigned int> bbIdsInThis;
|
| for (unsigned int i = 0; i < bbs_.size(); i++) {
|
| bbIdsInThis.push_back(bbs_[i]->getID());
|
| bbIdToThisBBIndex[bbs_[i]->getID()] = i;
|
| for (unsigned int j = 0; j < size_; j++) {
|
| if (bbs_[i]->getID() == other.bbs_[j]->getID()) {
|
| bbIdToOtherBBIndex[bbs_[i]->getID()] = j;
|
| break;
|
| }
|
| }
|
| }
|
|
|
|
|
| Match cmMatch;
|
| Molecule<Vector3> cmA, cmB;
|
|
|
| for (unsigned int bbId : bbIdsInThis) {
|
| Vector3 a = trans_[bbIdToThisBBIndex[bbId]] * bbs_[bbIdToThisBBIndex[bbId]]->cm_;
|
| cmA.add(a);
|
| Vector3 b = other.trans_[bbIdToOtherBBIndex[bbId]] * other.bbs_[bbIdToOtherBBIndex[bbId]]->cm_;
|
| cmB.add(b);
|
| }
|
| for (unsigned int i = 0; i < size_; i++) {
|
| cmMatch.add(i, i);
|
| }
|
| cmMatch.calculateBestFit(cmA, cmB);
|
| float foundRMSD = cmMatch.rmsd();
|
|
|
| int MAX_ITER = 10;
|
|
|
| for (int i = 0; i < MAX_ITER; i++) {
|
|
|
| bool somethingChanged = false;
|
| for (std::vector<unsigned int> identGroup : presentIdentGroups) {
|
| for (unsigned int i = 0; i < identGroup.size(); i++) {
|
| for (unsigned int j = i + 1; j < identGroup.size(); j++) {
|
| Match cmMatch;
|
| Molecule<Vector3> cmA, cmB;
|
|
|
| for (unsigned int bbId : bbIdsInThis) {
|
| Vector3 a;
|
| if (bbId == identGroup[i]) {
|
| unsigned int replacedbbId = identGroup[j];
|
| a = trans_[bbIdToThisBBIndex[replacedbbId]] * bbs_[bbIdToThisBBIndex[replacedbbId]]->cm_;
|
| } else if (bbId == identGroup[j]) {
|
| unsigned int replacedbbId = identGroup[i];
|
| a = trans_[bbIdToThisBBIndex[replacedbbId]] * bbs_[bbIdToThisBBIndex[replacedbbId]]->cm_;
|
| } else {
|
| a = trans_[bbIdToThisBBIndex[bbId]] * bbs_[bbIdToThisBBIndex[bbId]]->cm_;
|
| }
|
| cmA.add(a);
|
| Vector3 b =
|
| other.trans_[bbIdToOtherBBIndex[bbId]] * other.bbs_[bbIdToOtherBBIndex[bbId]]->cm_;
|
| cmB.add(b);
|
| }
|
| for (unsigned int i = 0; i < size_; i++) {
|
| cmMatch.add(i, i);
|
| }
|
| cmMatch.calculateBestFit(cmA, cmB);
|
| float newRMSD = cmMatch.rmsd();
|
| if (newRMSD < foundRMSD) {
|
|
|
|
|
|
|
| unsigned int tmp = bbIdToThisBBIndex[identGroup[i]];
|
| bbIdToThisBBIndex[identGroup[i]] = bbIdToThisBBIndex[identGroup[j]];
|
| bbIdToThisBBIndex[identGroup[j]] = tmp;
|
| foundRMSD = newRMSD;
|
| somethingChanged = true;
|
| } else {
|
|
|
|
|
| }
|
| }
|
| }
|
| }
|
| if (!somethingChanged)
|
| break;
|
| if (i == MAX_ITER - 1)
|
| std::cout << "stops clustering because max iteration " << bitIDS_ << std::endl;
|
| }
|
|
|
| if (foundRMSD > 1.5)
|
| return foundRMSD;
|
|
|
| Match match;
|
| Molecule<Atom> A, B;
|
| for (unsigned int bbId : bbIdsInThis) {
|
| Molecule<Atom> molA = bbs_[bbIdToThisBBIndex[bbId]]->caAtoms_;
|
| molA.rigidTrans(trans_[bbIdToThisBBIndex[bbId]]);
|
| A.concat(molA);
|
|
|
| Molecule<Atom> molB = other.bbs_[bbIdToOtherBBIndex[bbId]]->caAtoms_;
|
| molB.rigidTrans(other.trans_[bbIdToOtherBBIndex[bbId]]);
|
| B.concat(molB);
|
| }
|
|
|
| for (unsigned int i = 0; i < A.size(); i++) {
|
| match.add(i, i);
|
| }
|
| match.calculateBestFit(A, B);
|
| return match.rmsd();
|
| }
|
|
|
| double SuperBB::calcRmsd(const SuperBB &other) const {
|
| std::vector<std::vector<unsigned int>> emptyIdentGroups;
|
| return calcRmsd(other, emptyIdentGroups);
|
| }
|
|
|
| void SuperBB::fullReport(std::ostream &s) {
|
| HierarchicalFold::countResults_++;
|
| RigidTrans3 tr;
|
| s << "size_ " << size_ << " transScore_ " << transScore_ << " multPen_ " << multPen_ << " singlePen_ "
|
| << singlePen_ << " diffPen " << singlePen_ - multPen_ << " backBonePen_ " << backBonePen_ << " maxPen_ "
|
| << maxPen_ << " restraintsRatio_ " << restraintsRatio_ << " weightedTransScore " << weightedTransScore_;
|
| s << " [";
|
| int i = 0;
|
| for (auto it = bbs_.begin(); it != bbs_.end(); it++, i++) {
|
| s << (*it)->id_ << "(" << trans_[i] << ")";
|
| if (it + 1 != bbs_.end())
|
| s << ",";
|
| }
|
| s << "] " << tr;
|
| FoldStep::outputFoldSteps(s, foldSteps_);
|
| s << std::endl;
|
| }
|
|
|
| std::ostream &operator<<(std::ostream &s, const SuperBB &sbb) {
|
|
|
| s << "size_ " << sbb.size_ << " backBonePen_ " << sbb.backBonePen_ << " restraintsRatio_ " << sbb.restraintsRatio_;
|
|
|
|
|
| s << " [";
|
| int i = 0;
|
| for (auto it = sbb.bbs_.begin(); it != sbb.bbs_.end(); it++, i++) {
|
| s << (*it)->getID() << "(" << sbb.trans_[i] << ")";
|
| if (it + 1 != sbb.bbs_.end())
|
| s << ",";
|
| }
|
| s << "] ";
|
|
|
| FoldStep::outputFoldSteps(s, sbb.foldSteps_);
|
| s << std::endl;
|
| return s;
|
| }
|
|
|
| TransIterator2::TransIterator2(const SuperBB &bb1, const SuperBB &bb2, int pbb1, int pbb2)
|
| : bb1_(bb1), bb2_(bb2), index1_(0), index2_(0) {
|
|
|
| std::shared_ptr<const BB> b1 = bb1_.bbs_[index1_];
|
| while (b1->getID() != pbb1) {
|
| b1 = bb1_.bbs_[++index1_];
|
| }
|
| std::shared_ptr<const BB> b2 = bb2_.bbs_[index2_];
|
| while (b2->getID() != pbb2) {
|
| b2 = bb2_.bbs_[++index2_];
|
| }
|
|
|
| mediatorTrans_ = !bb2_.trans_[index2_];
|
| trans_ = &(b1->getTransformations(b2->getID()));
|
| it_ = trans_->begin();
|
| if (!isAtEnd()) {
|
| generateTransformation();
|
| }
|
| }
|
|
|
