From: Max Date: Tue, 27 Sep 2022 12:22:38 +0000 (+0300) Subject: Первая говноверсия говнопроливновых говноспиралей. X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=947262102c67fb6a85aa50925bf0b7f638aa8f8a;p=alexxy%2Fgromacs-dssp.git Первая говноверсия говнопроливновых говноспиралей. Работать, конечно же, не будет. --- diff --git a/src/dssp.cpp b/src/dssp.cpp index 2818c00..00587fb 100644 --- a/src/dssp.cpp +++ b/src/dssp.cpp @@ -103,6 +103,8 @@ void Dssp::initOptions(IOptionsContainer* options, TrajectoryAnalysisSettings* s BooleanOption("addh").store(&initParams.addHydrogens).defaultValue(false).description("Add hydrogens or not")); options->addOption( BooleanOption("p").store(&initParams.PPHelices).defaultValue(true).description("Prefer Pi Helices")); + options->addOption( + IntegerOption("s").store(&initParams.pp_stretch).defaultValue(3).description("Stretch for your pp. 2 or 3")); options->addOption( BooleanOption("v").store(&initParams.verbose).defaultValue(false).description(">:(")); settings->setHelpText(desc); diff --git a/src/dssptools.cpp b/src/dssptools.cpp index bd3c2da..def60a2 100644 --- a/src/dssptools.cpp +++ b/src/dssptools.cpp @@ -82,19 +82,25 @@ std::size_t ResInfo::getIndex(backboneAtomTypes atomTypeName) const{ secondaryStructures::secondaryStructures(){ } -void secondaryStructures::initiateSearch(const std::vector &ResInfoMatrix, const bool PiHelicesPreferencez){ +void secondaryStructures::initiateSearch(const std::vector &ResInfoMatrix, const bool PiHelicesPreferencez, const int _pp_stretch){ SecondaryStructuresStatusMap.resize(0); SecondaryStructuresStringLine.resize(0); std::vector temp; temp.resize(0), PiHelixPreference = PiHelicesPreferencez; ResInfoMap = &ResInfoMatrix; + pp_stretch = _pp_stretch; SecondaryStructuresStatusMap.resize(ResInfoMatrix.size()); SecondaryStructuresStringLine.resize(ResInfoMatrix.size(), '~'); } -void secondaryStructures::secondaryStructuresData::setStatus(const secondaryStructureTypes secondaryStructureTypeName){ - SecondaryStructuresStatusArray[static_cast(secondaryStructureTypeName)] = true; - SecondaryStructuresStatus = secondaryStructureTypeName; +void secondaryStructures::secondaryStructuresData::setStatus(const secondaryStructureTypes secondaryStructureTypeName, bool status){ + SecondaryStructuresStatusArray[static_cast(secondaryStructureTypeName)] = status; + if (status){ + SecondaryStructuresStatus = secondaryStructureTypeName; + } + else{ + SecondaryStructuresStatus = secondaryStructureTypes::Loop; + } } void secondaryStructures::secondaryStructuresData::setStatus(const HelixPositions helixPosition, const turnsTypes turn){ @@ -291,10 +297,6 @@ void secondaryStructures::analyzeBridgesAndStrandsPatterns(){ for(const bridgeTypes &bridgeType : {bridgeTypes::ParallelBridge, bridgeTypes::AntiParallelBridge}){ if (SecondaryStructuresStatusMap[i].hasBridges(bridgeType) && SecondaryStructuresStatusMap[i + j].hasBridges(bridgeType) ){ std::size_t i_partner{SecondaryStructuresStatusMap[i].getBridgePartnerIndex(bridgeType)}, j_partner{SecondaryStructuresStatusMap[i + j].getBridgePartnerIndex(bridgeType)}, second_strand{}; -// std::cout << "i = " << i << std::endl; -// std::cout << "j = " << j << std::endl; -// std::cout << "i_partner = " << i_partner << std::endl; -// std::cout << "j_partner = " << j_partner << std::endl; if ( abs(i_partner - j_partner) < 6){ if (i_partner < j_partner){ second_strand = i_partner; @@ -303,18 +305,22 @@ void secondaryStructures::analyzeBridgesAndStrandsPatterns(){ second_strand = j_partner; } for(int k{abs(i_partner - j_partner)}; k >= 0; --k){ + if (SecondaryStructuresStatusMap[second_strand + k].getStatus(secondaryStructureTypes::Bridge)){ + SecondaryStructuresStatusMap[second_strand + k].setStatus(secondaryStructureTypes::Bridge, false); + } + SecondaryStructuresStatusMap[second_strand + k].setStatus(secondaryStructureTypes::Strand); - std::cout << second_strand + k << " is strand" << std::endl; } is_Estrand = true; } - } } if (is_Estrand){ for(std::size_t k{0}; k <= j; ++k){ + if (SecondaryStructuresStatusMap[i + k].getStatus(secondaryStructureTypes::Bridge)){ + SecondaryStructuresStatusMap[i + k].setStatus(secondaryStructureTypes::Bridge, false); + } SecondaryStructuresStatusMap[i + k].setStatus(secondaryStructureTypes::Strand); - std::cout << i + k << " is strand" << std::endl; } break; } @@ -389,11 +395,8 @@ void secondaryStructures::analyzeBridgesAndStrandsPatterns(){ void secondaryStructures::analyzeTurnsAndHelicesPatterns(){ for(const turnsTypes &i : { turnsTypes::Turn_4, turnsTypes::Turn_3, turnsTypes::Turn_5 }){ std::size_t stride {static_cast(i) + 3}; -// std::cout << "Testing Helix_" << stride << std::endl; for(std::size_t j {0}; j + stride < SecondaryStructuresStatusMap.size(); ++j){ -// std::cout << "Testing " << j << " and " << j + stride << std::endl; if ( hasHBondBetween(j + stride, j) && NoChainBreaksBetween(j, j + stride) ){ -// std::cout << j << " and " << j + stride << " has hbond!" << std::endl; SecondaryStructuresStatusMap[j + stride].setStatus(HelixPositions::End, i); for (std::size_t k {1}; k < stride; ++k){ @@ -465,14 +468,11 @@ void secondaryStructures::analyzeTurnsAndHelicesPatterns(){ } -void secondaryStructures::analyzePPHelicesPatterns(){} - std::string secondaryStructures::patternSearch(){ analyzeBridgesAndStrandsPatterns(); analyzeTurnsAndHelicesPatterns(); -// analyzePPHelicesPatterns(); // for(std::size_t i {0}; i < ResInfoMap->size(); ++i){ // std::cout << (*ResInfoMap)[i].info->nr << " " << *((*ResInfoMap)[i].info->name) << std::endl; @@ -738,6 +738,148 @@ void DsspTool::calculateBends(const t_trxframe &fr, const t_pbc *pbc) } } +float DsspTool::CalculateDihedralAngle(const int &A, const int &B, const int &C, const int &D, const t_trxframe &fr, const t_pbc *pbc){ + float result{360}, u{}, v{}; + gmx::RVec vBA{}, vDC{}, vBD{}, p{}, x{}, y{}; + pbc_dx(pbc, fr.x[A], fr.x[B], vBA.as_vec()); + pbc_dx(pbc, fr.x[C], fr.x[D], vDC.as_vec()); + pbc_dx(pbc, fr.x[D], fr.x[B], vBD.as_vec()); + + for(std::size_t i{XX}, j{j + 1}, k{i + 2}; i <= ZZ; ++i, ++j, ++k){ + if (j > 2){ + j -= 3; + } + if (k > 2){ + k -= 3; + } + p[i] = (vBD[j] * vBA[k]) - (vBD[k] * vBA[j]); + x[i] = (vBD[j] * vDC[k]) - (vBD[k] * vDC[j]); + + } + for(std::size_t i{XX}, j{j + 1}, k{i + 2}; i <= ZZ; ++i, ++j, ++k){ + if (j > 2){ + j -= 3; + } + if (k > 2){ + k -= 3; + } + y[i] = (vBD[j] * x[k]) - (vBD[k] * x[j]); + } + + u = (x[XX] * x[XX]) + (x[YY] * x[YY]) + (x[ZZ] * x[ZZ]); + v = (y[XX] * y[XX]) + (y[YY] * y[YY]) + (y[ZZ] * y[ZZ]); + + if (u > 0 and v > 0){ + u = ((p[XX] * x[XX]) + (p[YY] * x[YY]) + (p[ZZ] * x[ZZ])) / std::sqrt(u); + v = ((p[XX] * y[XX]) + (p[YY] * y[YY]) + (p[ZZ] * y[ZZ])) / std::sqrt(v); + if (u != 0 or v != 0){ + result = std::atan2(v, u) * gmx::c_rad2Deg; + } + } + return result; +} + +void DsspTool::calculateDihedrals(const t_trxframe &fr, const t_pbc *pbc){ + const float epsilon = 29; + const float phi_min = -75 - epsilon; + const float phi_max = -75 + epsilon; + const float psi_min = 145 - epsilon; + const float psi_max = 145 + epsilon; + std::vector phi, psi; + phi.resize(0); + psi.resize(0); + phi.resize(IndexMap.size(), 360); + psi.resize(IndexMap.size(), 360); + + for (std::size_t i = 1; i + 1 < IndexMap.size(); ++i){ // TODO add index verifictaion (check if those atom indexes exist) + phi[i] = CalculateDihedralAngle(static_cast(IndexMap[i - 1].getIndex(backboneAtomTypes::AtomC)), + static_cast(IndexMap[i].getIndex(backboneAtomTypes::AtomN)), + static_cast(IndexMap[i].getIndex(backboneAtomTypes::AtomCA)), + static_cast(IndexMap[i].getIndex(backboneAtomTypes::AtomC)), + fr, + pbc); + psi[i] = CalculateDihedralAngle(static_cast(IndexMap[i].getIndex(backboneAtomTypes::AtomN)), + static_cast(IndexMap[i].getIndex(backboneAtomTypes::AtomCA)), + static_cast(IndexMap[i].getIndex(backboneAtomTypes::AtomC)), + static_cast(IndexMap[i + 1].getIndex(backboneAtomTypes::AtomN)), + fr, + pbc); + } + + for (std::size_t i = 1; i + 3 < IndexMap.size(); ++i){ + switch (initParams.pp_stretch){ + case 2: { + if (phi_min > phi[i] or phi[i] > phi_max or + phi_min > phi[i + 1] or phi[i + 1]> phi_max){ + continue; + } + + if (psi_min > psi[i] or psi[i] > psi_max or + psi_min > psi[i + 1] or psi[i + 1] > psi_max){ + continue; + } + + switch (PatternSearch.SecondaryStructuresStatusMap[i].getStatus(turnsTypes::Turn_PP)){ + case HelixPositions::None: + PatternSearch.SecondaryStructuresStatusMap[i].setStatus(HelixPositions::Start, turnsTypes::Turn_PP); + break; + + case HelixPositions::End: + PatternSearch.SecondaryStructuresStatusMap[i].setStatus(HelixPositions::Start_AND_End, turnsTypes::Turn_PP); + break; + + default: + break; + } + + PatternSearch.SecondaryStructuresStatusMap[i + 1].setStatus(HelixPositions::End, turnsTypes::Turn_PP); + /* Пропустил проверку того, что заменяемая ак - петля */ + PatternSearch.SecondaryStructuresStatusMap[i].setStatus(secondaryStructureTypes::Helix_PP); + PatternSearch.SecondaryStructuresStatusMap[i + 1].setStatus(secondaryStructureTypes::Helix_PP); + break; + } + case 3:{ + if (phi_min > phi[i] or phi[i] > phi_max or + phi_min > phi[i + 1] or phi[i + 1]> phi_max or + phi_min > phi[i + 2] or phi[i + 2]> phi_max){ + continue; + } + + if (psi_min > psi[i] or psi[i] > psi_max or + psi_min > psi[i + 1] or psi[i + 1] > psi_max or + psi_min > psi[i + 2] or psi[i + 2] > psi_max){ + continue; + } + + switch (PatternSearch.SecondaryStructuresStatusMap[i].getStatus(turnsTypes::Turn_PP)){ + case HelixPositions::None: + PatternSearch.SecondaryStructuresStatusMap[i].setStatus(HelixPositions::Start, turnsTypes::Turn_PP); + break; + + case HelixPositions::End: + PatternSearch.SecondaryStructuresStatusMap[i].setStatus(HelixPositions::Start_AND_End, turnsTypes::Turn_PP); + break; + + default: + break; + } + + PatternSearch.SecondaryStructuresStatusMap[i + 1].setStatus(HelixPositions::Middle, turnsTypes::Turn_PP); + PatternSearch.SecondaryStructuresStatusMap[i + 2].setStatus(HelixPositions::End, turnsTypes::Turn_PP); + /* Пропустил проверку того, что заменяемая ак - петля */ + PatternSearch.SecondaryStructuresStatusMap[i].setStatus(secondaryStructureTypes::Helix_PP); + PatternSearch.SecondaryStructuresStatusMap[i + 1].setStatus(secondaryStructureTypes::Helix_PP); + PatternSearch.SecondaryStructuresStatusMap[i + 2].setStatus(secondaryStructureTypes::Helix_PP); + + break; + } + default: + throw std::runtime_error("Unsupported stretch length"); + } + } + +} + void DsspTool::calculateHBondEnergy(ResInfo& Donor, ResInfo& Acceptor, const t_trxframe& fr, @@ -1024,7 +1166,7 @@ void DsspTool::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc) // std::cout << IndexMap[i].info->nr << " " << *(IndexMap[i].info->name) << std::endl; // } - PatternSearch.initiateSearch(IndexMap, initParams.PPHelices); + PatternSearch.initiateSearch(IndexMap, initParams.PPHelices, initParams.pp_stretch); calculateBends(fr, pbc); Storage.storageData(frnr, PatternSearch.patternSearch()); diff --git a/src/dssptools.h b/src/dssptools.h index 048d153..940dd42 100644 --- a/src/dssptools.h +++ b/src/dssptools.h @@ -77,6 +77,7 @@ struct initParameters { real cutoff_; // = 4.0; ??? NBSearchMethod NBS; bool verbose, PPHelices, addHydrogens; + int pp_stretch; }; enum class backboneAtomTypes : std::size_t @@ -122,11 +123,11 @@ enum class secondaryStructureTypes : std::size_t { // TODO Break, // = Bend, // S Turn, // T - Helis_PPII, // P + Helix_PP, // P Helix_5, // I Helix_3, // G - Bridge, // B Strand, // E + Bridge, // B Helix_4, // H Count @@ -136,6 +137,7 @@ enum class turnsTypes : std::size_t { Turn_3 = 0, Turn_4, Turn_5, + Turn_PP, Count }; @@ -158,13 +160,13 @@ enum class bridgeTypes : std::size_t { class secondaryStructures{ // PatterSearch Wrapper public: secondaryStructures(); - void initiateSearch(const std::vector &ResInfoMatrix, const bool PiHelicesPreferencez = true); + void initiateSearch(const std::vector &ResInfoMatrix, const bool PiHelicesPreferencez, const int _pp_stretch); std::string patternSearch(); ~secondaryStructures(); class secondaryStructuresData{ // PatternSearch Tool public: - void setStatus(const secondaryStructureTypes secondaryStructureTypeName); + void setStatus(const secondaryStructureTypes secondaryStructureTypeName, bool status = true); void setStatus(const HelixPositions helixPosition, const turnsTypes turn); bool getStatus(const secondaryStructureTypes secondaryStructureTypeName) const, isBreakPartnerWith(const secondaryStructuresData *partner) const; HelixPositions getStatus(const turnsTypes turn) const; @@ -173,13 +175,14 @@ public: bool hasBridges() const, hasBridges(bridgeTypes bridgeType) const, isBridgePartnerWith(secondaryStructuresData *bridgePartner, bridgeTypes bridgeType) const; std::size_t getBridgePartnerIndex(bridgeTypes bridgeType) const; secondaryStructuresData getBridgePartner(bridgeTypes bridgeType) const; + void set_phi(const ResInfo resi), set_psi(const ResInfo resi); private: std::array(secondaryStructureTypes::Count)> SecondaryStructuresStatusArray{ true, 0, 0, 0, 0, 0, 0 }; secondaryStructuresData *breakPartners[2]{nullptr, nullptr}; secondaryStructuresData *bridgePartners[2]{nullptr, nullptr}; std::size_t bridgePartnersIndexes[2]{}; secondaryStructureTypes SecondaryStructuresStatus {secondaryStructureTypes::Loop}; - std::array(turnsTypes::Count)> TurnsStatusArray {HelixPositions::None, HelixPositions::None, HelixPositions::None}; + std::array(turnsTypes::Count)> TurnsStatusArray {HelixPositions::None, HelixPositions::None, HelixPositions::None, HelixPositions::None}; }; std::vector SecondaryStructuresStatusMap; @@ -188,13 +191,14 @@ private: const std::vector *ResInfoMap; const gmx::EnumerationArray secondaryStructureTypeNames = { - { '~', '=', 'S', 'T', 'P', 'I', 'G', 'B', 'E', 'H'} + { '~', '=', 'S', 'T', 'P', 'I', 'G', 'E', 'B', 'H'} }; std::string SecondaryStructuresStringLine; bool hasHBondBetween(std::size_t resi1, std::size_t resi2) const; bool NoChainBreaksBetween(std::size_t ResiStart, std::size_t ResiEnd) const, isLoop(const std::size_t resiIndex) const, PiHelixPreference; + int8_t pp_stretch; bridgeTypes calculateBridge(std::size_t i, std::size_t j) const; void analyzeBridgesAndStrandsPatterns(), analyzeTurnsAndHelicesPatterns(), analyzePPHelicesPatterns(); @@ -239,7 +243,7 @@ private: initParameters initParams; AtomsDataPtr atoms_; std::string filename_; - void calculateBends(const t_trxframe& fr, const t_pbc* pbc); + void calculateBends(const t_trxframe& fr, const t_pbc* pbc), calculateDihedrals(const t_trxframe& fr, const t_pbc* pbc); std::vector ResiNames; std::vector AtomResi; std::vector IndexMap; @@ -253,6 +257,7 @@ private: const t_pbc* pbc); float CalculateAtomicDistances(const int &A, const int &B, const t_trxframe &fr, const t_pbc *pbc); float CalculateAtomicDistances(const rvec &A, const int &B, const t_trxframe &fr, const t_pbc *pbc); + float CalculateDihedralAngle(const int &A, const int &B, const int &C, const int &D, const t_trxframe &fr, const t_pbc *pbc); class DsspStorage {