| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | /* BEGIN software license | ||
| 2 | * | ||
| 3 | * MsXpertSuite - mass spectrometry software suite | ||
| 4 | * ----------------------------------------------- | ||
| 5 | * Copyright(C) 2009,...,2018 Filippo Rusconi | ||
| 6 | * | ||
| 7 | * http://www.msxpertsuite.org | ||
| 8 | * | ||
| 9 | * This file is part of the MsXpertSuite project. | ||
| 10 | * | ||
| 11 | * The MsXpertSuite project is the successor of the massXpert project. This | ||
| 12 | * project now includes various independent modules: | ||
| 13 | * | ||
| 14 | * - massXpert, model polymer chemistries and simulate mass spectrometric data; | ||
| 15 | * - mineXpert, a powerful TIC chromatogram/mass spectrum viewer/miner; | ||
| 16 | * | ||
| 17 | * This program is free software: you can redistribute it and/or modify | ||
| 18 | * it under the terms of the GNU General Public License as published by | ||
| 19 | * the Free Software Foundation, either version 3 of the License, or | ||
| 20 | * (at your option) any later version. | ||
| 21 | * | ||
| 22 | * This program is distributed in the hope that it will be useful, | ||
| 23 | * but WITHOUT ANY WARRANTY; without even the implied warranty of | ||
| 24 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | ||
| 25 | * GNU General Public License for more details. | ||
| 26 | * | ||
| 27 | * You should have received a copy of the GNU General Public License | ||
| 28 | * along with this program. If not, see <http://www.gnu.org/licenses/>. | ||
| 29 | * | ||
| 30 | * END software license | ||
| 31 | */ | ||
| 32 | |||
| 33 | |||
| 34 | /////////////////////// Std includes | ||
| 35 | #include <cmath> | ||
| 36 | |||
| 37 | |||
| 38 | /////////////////////// Local includes | ||
| 39 | #include "MsXpS/libXpertMassCore/PkaPhPi.hpp" | ||
| 40 | |||
| 41 | |||
| 42 | namespace MsXpS | ||
| 43 | { | ||
| 44 | namespace libXpertMassCore | ||
| 45 | { | ||
| 46 | |||
| 47 | |||
| 48 | /*! | ||
| 49 | \class MsXpS::libXpertMassCore::PkaPhPi | ||
| 50 | \inmodule libXpertMassCore | ||
| 51 | \ingroup PolChemDefBuildingdBlocks | ||
| 52 | \inheaderfile PkaPhPi.hpp | ||
| 53 | |||
| 54 | \brief The PkaPhPi class provides a model for specifying the | ||
| 55 | acido-basic properties of a chemical entity. | ||
| 56 | */ | ||
| 57 | |||
| 58 | |||
| 59 | /*! | ||
| 60 | \variable MsXpS::libXpertMassCore::PkaPhPi::mcsp_polChemDef | ||
| 61 | |||
| 62 | \brief The polymer chemistry definition. | ||
| 63 | */ | ||
| 64 | |||
| 65 | /*! | ||
| 66 | \variable MsXpS::libXpertMassCore::PkaPhPi::m_monomers | ||
| 67 | |||
| 68 | \brief The container of \l Monomer instances as read from the pka_ph_pi.xml | ||
| 69 | file. | ||
| 70 | */ | ||
| 71 | |||
| 72 | /*! | ||
| 73 | \variable MsXpS::libXpertMassCore::PkaPhPi::m_modifs | ||
| 74 | |||
| 75 | \brief The container of \l Modif instances as read from the pka_ph_pi.xml file. | ||
| 76 | */ | ||
| 77 | |||
| 78 | /*! | ||
| 79 | \variable MsXpS::libXpertMassCore::PkaPhPi::m_ph | ||
| 80 | |||
| 81 | \brief The pH of the environment. | ||
| 82 | |||
| 83 | This pH value is required to compute the number of charges of a given chemical | ||
| 84 | entity (a \l Polymer) sequence, for example. | ||
| 85 | */ | ||
| 86 | |||
| 87 | /*! | ||
| 88 | \variable MsXpS::libXpertMassCore::PkaPhPi::m_pi | ||
| 89 | |||
| 90 | \brief The pI of the chemical entity. | ||
| 91 | */ | ||
| 92 | |||
| 93 | /*! | ||
| 94 | \variable MsXpS::libXpertMassCore::PkaPhPi::mcsp_polymer | ||
| 95 | |||
| 96 | \brief The polymer of which the acidobasic properties are computed. | ||
| 97 | */ | ||
| 98 | |||
| 99 | /*! | ||
| 100 | \variable MsXpS::libXpertMassCore::PkaPhPi::m_calcOptions | ||
| 101 | |||
| 102 | \brief The \l CalcOptions that configure the way the computations are to be | ||
| 103 | carried out. | ||
| 104 | */ | ||
| 105 | |||
| 106 | /*! | ||
| 107 | \variable MsXpS::libXpertMassCore::PkaPhPi::m_positiveCharges | ||
| 108 | |||
| 109 | \brief The count of positive charges. | ||
| 110 | */ | ||
| 111 | |||
| 112 | /*! | ||
| 113 | \variable MsXpS::libXpertMassCore::PkaPhPi::m_negativeCharges | ||
| 114 | |||
| 115 | \brief The count of negative charges. | ||
| 116 | */ | ||
| 117 | |||
| 118 | /*! | ||
| 119 | \brief Constructs a PkaPhPi instance with a number of parameters. | ||
| 120 | |||
| 121 | \list | ||
| 122 | |||
| 123 | \li \a pol_chem_def_csp: The polymer chemistry definition | ||
| 124 | |||
| 125 | \li \a polymer_cqsp: The polymer within the context of which the calculations are | ||
| 126 | performed. | ||
| 127 | |||
| 128 | \li \a calc_options: The options driving the calculations. | ||
| 129 | |||
| 130 | \endlist | ||
| 131 | */ | ||
| 132 | ✗ | PkaPhPi::PkaPhPi(PolChemDefCstSPtr pol_chem_def_csp, | |
| 133 | PolymerCstQSPtr polymer_cqsp, | ||
| 134 | ✗ | const CalcOptions &calc_options) | |
| 135 | ✗ | : mcsp_polChemDef(pol_chem_def_csp), | |
| 136 | ✗ | mcsp_polymer(polymer_cqsp), | |
| 137 | ✗ | m_calcOptions(calc_options) | |
| 138 | { | ||
| 139 | ✗ | if(mcsp_polChemDef == nullptr || mcsp_polChemDef.get() == nullptr) | |
| 140 | ✗ | qCritical() | |
| 141 | ✗ | << "Creating PkaPhPi instance without polymer chemistry definition."; | |
| 142 | |||
| 143 | ✗ | if(polymer_cqsp == nullptr) | |
| 144 | ✗ | qFatalStream() << "Programming error. Pointer cannot be nullptr."; | |
| 145 | ✗ | } | |
| 146 | |||
| 147 | /*! | ||
| 148 | \brief Destructs this PkaPhPi instance. | ||
| 149 | */ | ||
| 150 | ✗ | PkaPhPi::~PkaPhPi() | |
| 151 | { | ||
| 152 | ✗ | } | |
| 153 | |||
| 154 | /*! | ||
| 155 | \brief Sets the PolChemDef to \a pol_chem_def_csp. | ||
| 156 | */ | ||
| 157 | void | ||
| 158 | ✗ | PkaPhPi::setPolChemDefCstSPtr(PolChemDefCstSPtr pol_chem_def_csp) | |
| 159 | { | ||
| 160 | ✗ | mcsp_polChemDef = pol_chem_def_csp; | |
| 161 | |||
| 162 | ✗ | if(mcsp_polChemDef == nullptr || mcsp_polChemDef.get() == nullptr) | |
| 163 | ✗ | qCritical() << "Setting an Creating PkaPhPi instance without polymer " | |
| 164 | ✗ | "chemistry definition."; | |
| 165 | ✗ | } | |
| 166 | |||
| 167 | /*! | ||
| 168 | \brief Returns the PolChemDef. | ||
| 169 | */ | ||
| 170 | PolChemDefCstSPtr | ||
| 171 | ✗ | PkaPhPi::getPolChemDefCstSPtr() const | |
| 172 | { | ||
| 173 | ✗ | return mcsp_polChemDef; | |
| 174 | } | ||
| 175 | |||
| 176 | /*! | ||
| 177 | \brief Moves the Monomer instances from \a monomers to the member container. | ||
| 178 | */ | ||
| 179 | void | ||
| 180 | ✗ | PkaPhPi::setMonomers(std::vector<MonomerCstSPtr> &&monomers) | |
| 181 | { | ||
| 182 | ✗ | m_monomers.clear(); | |
| 183 | ✗ | m_monomers = std::move(monomers); | |
| 184 | ✗ | } | |
| 185 | |||
| 186 | /*! | ||
| 187 | \brief Moves the Monomer instances from \a monomers to the member container. | ||
| 188 | */ | ||
| 189 | void | ||
| 190 | ✗ | PkaPhPi::setMonomers(std::vector<MonomerSPtr> &&monomers) | |
| 191 | { | ||
| 192 | // We want the pointers to be to const Monomer. | ||
| 193 | |||
| 194 | ✗ | m_monomers.clear(); | |
| 195 | ✗ | m_monomers.reserve(monomers.size()); | |
| 196 | |||
| 197 | ✗ | std::transform(std::make_move_iterator(monomers.begin()), | |
| 198 | std::make_move_iterator(monomers.end()), | ||
| 199 | ✗ | std::back_inserter(m_monomers), | |
| 200 | ✗ | [](std::shared_ptr<Monomer> &&monomer_sp) { | |
| 201 | ✗ | return std::const_pointer_cast<const Monomer>( | |
| 202 | ✗ | std::move(monomer_sp)); | |
| 203 | }); | ||
| 204 | |||
| 205 | ✗ | monomers.clear(); | |
| 206 | ✗ | } | |
| 207 | |||
| 208 | /*! | ||
| 209 | \brief Copies the Monomer instances in \a monomers to the member container. | ||
| 210 | |||
| 211 | This is a shallow copy with only the pointers being copied. | ||
| 212 | */ | ||
| 213 | void | ||
| 214 | ✗ | PkaPhPi::setMonomers(std::vector<MonomerCstSPtr> &monomers) | |
| 215 | { | ||
| 216 | ✗ | m_monomers = monomers; | |
| 217 | ✗ | } | |
| 218 | |||
| 219 | |||
| 220 | /*! | ||
| 221 | \brief Returns a const reference to the container of Monomer instances. | ||
| 222 | */ | ||
| 223 | const std::vector<MonomerCstSPtr> & | ||
| 224 | ✗ | PkaPhPi::getMonomersCstRef() const | |
| 225 | { | ||
| 226 | ✗ | return m_monomers; | |
| 227 | } | ||
| 228 | |||
| 229 | /*! | ||
| 230 | \brief Moves the Modif instances from \a modifs to the member container. | ||
| 231 | */ | ||
| 232 | void | ||
| 233 | ✗ | PkaPhPi::setModifs(std::vector<ModifCstSPtr> &&modifs) | |
| 234 | { | ||
| 235 | ✗ | m_modifs.clear(); | |
| 236 | ✗ | m_modifs = std::move(modifs); | |
| 237 | ✗ | modifs.clear(); | |
| 238 | ✗ | } | |
| 239 | |||
| 240 | /*! | ||
| 241 | \brief Moves the Modif instances from \a modifs to the member container. | ||
| 242 | */ | ||
| 243 | void | ||
| 244 | ✗ | PkaPhPi::setModifs(std::vector<ModifSPtr> &&modifs) | |
| 245 | { | ||
| 246 | // We want the pointers to be to const Modif. | ||
| 247 | |||
| 248 | ✗ | m_modifs.clear(); | |
| 249 | ✗ | m_modifs.reserve(modifs.size()); | |
| 250 | |||
| 251 | ✗ | std::transform(std::make_move_iterator(modifs.begin()), | |
| 252 | std::make_move_iterator(modifs.end()), | ||
| 253 | ✗ | std::back_inserter(m_modifs), | |
| 254 | ✗ | [](std::shared_ptr<Modif> &&modif_sp) { | |
| 255 | ✗ | return std::const_pointer_cast<const Modif>( | |
| 256 | ✗ | std::move(modif_sp)); | |
| 257 | }); | ||
| 258 | |||
| 259 | ✗ | modifs.clear(); | |
| 260 | ✗ | } | |
| 261 | |||
| 262 | /*! | ||
| 263 | \brief Copies the Modif instances in \a modifs to the member container. | ||
| 264 | |||
| 265 | This is a shallow copy with only the pointers being copied. | ||
| 266 | */ | ||
| 267 | void | ||
| 268 | ✗ | PkaPhPi::setModifs(std::vector<ModifCstSPtr> &modifs) | |
| 269 | { | ||
| 270 | ✗ | m_modifs = modifs; | |
| 271 | ✗ | } | |
| 272 | |||
| 273 | /*! | ||
| 274 | \brief Returns a const reference to the container of Modif instances. | ||
| 275 | */ | ||
| 276 | const std::vector<ModifCstSPtr> & | ||
| 277 | ✗ | PkaPhPi::getModifs() const | |
| 278 | { | ||
| 279 | ✗ | return m_modifs; | |
| 280 | } | ||
| 281 | |||
| 282 | |||
| 283 | /*! | ||
| 284 | \brief Sets the pH to \a ph. | ||
| 285 | */ | ||
| 286 | void | ||
| 287 | ✗ | PkaPhPi::setPh(double ph) | |
| 288 | { | ||
| 289 | ✗ | Q_ASSERT(ph > 0 && ph < 14); | |
| 290 | |||
| 291 | ✗ | m_ph = ph; | |
| 292 | ✗ | } | |
| 293 | |||
| 294 | |||
| 295 | /*! | ||
| 296 | \brief Returns the pH. | ||
| 297 | */ | ||
| 298 | double | ||
| 299 | ✗ | PkaPhPi::ph() | |
| 300 | { | ||
| 301 | ✗ | return m_ph; | |
| 302 | } | ||
| 303 | |||
| 304 | /*! | ||
| 305 | \brief Returns the pI. | ||
| 306 | */ | ||
| 307 | double | ||
| 308 | ✗ | PkaPhPi::pi() | |
| 309 | { | ||
| 310 | ✗ | return m_pi; | |
| 311 | } | ||
| 312 | |||
| 313 | /*! | ||
| 314 | \brief Returns the positive charges. | ||
| 315 | */ | ||
| 316 | double | ||
| 317 | ✗ | PkaPhPi::positiveCharges() | |
| 318 | { | ||
| 319 | ✗ | return m_positiveCharges; | |
| 320 | } | ||
| 321 | |||
| 322 | /*! | ||
| 323 | \brief Returns the negative charges. | ||
| 324 | */ | ||
| 325 | double | ||
| 326 | ✗ | PkaPhPi::negativeCharges() | |
| 327 | { | ||
| 328 | ✗ | return m_negativeCharges; | |
| 329 | } | ||
| 330 | |||
| 331 | /*! | ||
| 332 | \brief Sets the calculation options to \a calc_options. | ||
| 333 | */ | ||
| 334 | void | ||
| 335 | ✗ | PkaPhPi::setCalcOptions(const libXpertMassCore::CalcOptions &calc_options) | |
| 336 | { | ||
| 337 | ✗ | m_calcOptions.initialize(calc_options); | |
| 338 | ✗ | } | |
| 339 | |||
| 340 | /*! | ||
| 341 | \brief Calculates the charges (positive and negative). | ||
| 342 | |||
| 343 | The general scheme is : | ||
| 344 | |||
| 345 | Get the list of the iter_index_range of the different \l Polymer region | ||
| 346 | selections. For each first monomer and end monomer of a given region selection, | ||
| 347 | check if the the region is an oligomer or a residual chain (m_selectionType of | ||
| 348 | libXpertMassCore::CalcOptions); act accordingly. Also, check for each selection | ||
| 349 | region if it encompasses the polymer left/right end. If the left/right end | ||
| 350 | modifications are to be taken into account, act accordingly. | ||
| 351 | |||
| 352 | The positive and negative charges are stored in the member \l m_positiveCharges | ||
| 353 | and \l m_negativeCharges variables. | ||
| 354 | |||
| 355 | Returns the count of chemical groups that have been processed. | ||
| 356 | |||
| 357 | \sa calculatePi() | ||
| 358 | */ | ||
| 359 | int | ||
| 360 | ✗ | PkaPhPi::calculateCharges() | |
| 361 | { | ||
| 362 | ✗ | int processedChemicalGroups = 0; | |
| 363 | |||
| 364 | ✗ | m_positiveCharges = 0; | |
| 365 | ✗ | m_negativeCharges = 0; | |
| 366 | |||
| 367 | // We of course need monomers ! Instead, we may not need modifs. | ||
| 368 | ✗ | if(mcsp_polChemDef == nullptr || mcsp_polChemDef.get() == nullptr) | |
| 369 | return -1; | ||
| 370 | |||
| 371 | ✗ | std::size_t polymerSize = mcsp_polymer->size(); | |
| 372 | |||
| 373 | ✗ | const IndexRangeCollection &index_range_collection = | |
| 374 | ✗ | m_calcOptions.getIndexRangeCollectionCstRef(); | |
| 375 | |||
| 376 | ✗ | for(qsizetype iter = 0; iter < index_range_collection.size(); ++iter) | |
| 377 | { | ||
| 378 | ✗ | const IndexRange &iter_index_range = | |
| 379 | ✗ | index_range_collection.getRangeCstRefAt(iter); | |
| 380 | |||
| 381 | ✗ | qsizetype start_index = iter_index_range.m_start; | |
| 382 | ✗ | qsizetype stop_index = iter_index_range.m_stop; | |
| 383 | |||
| 384 | ✗ | bool is_left_most_sequence_range = | |
| 385 | ✗ | index_range_collection.isLeftMostIndexRange(iter_index_range); | |
| 386 | ✗ | bool is_right_most_sequence_range = | |
| 387 | ✗ | index_range_collection.isRightMostIndexRange(iter_index_range); | |
| 388 | |||
| 389 | ✗ | for(qsizetype jter = start_index; jter < stop_index + 1; ++jter) | |
| 390 | { | ||
| 391 | ✗ | MonomerSPtr monomer_csp = | |
| 392 | ✗ | mcsp_polymer->getSequenceCstRef().getMonomerCstSPtrAt(jter); | |
| 393 | |||
| 394 | // Find a monomer by the same code in our list of monomers | ||
| 395 | // that have been fed with chemical group data. Note that | ||
| 396 | // all the monomers in a given sequence must not | ||
| 397 | // necessarily have a counterpart in the local list of | ||
| 398 | // monoemers. For example, there might be cases in which a | ||
| 399 | // given monomer might not bring any charge whatsoever. | ||
| 400 | |||
| 401 | ✗ | QString code = monomer_csp->getCode(); | |
| 402 | |||
| 403 | ✗ | std::vector<MonomerCstSPtr>::const_iterator monomer_iterator_cst = | |
| 404 | ✗ | std::find_if(m_monomers.cbegin(), | |
| 405 | m_monomers.cend(), | ||
| 406 | ✗ | [code](const MonomerCstSPtr &monomer_csp) { | |
| 407 | ✗ | return monomer_csp->getCode() == code; | |
| 408 | }); | ||
| 409 | |||
| 410 | ✗ | if(monomer_iterator_cst == m_monomers.cend()) | |
| 411 | ✗ | continue; | |
| 412 | |||
| 413 | // A monomer can have multiple such "CHEMICAL_GROUP" | ||
| 414 | // properties. Indeed, for example for proteins, a monomer | ||
| 415 | // might have three such chemical groups(and thus three | ||
| 416 | // Prop objects): one for the alpha NH2, one for the | ||
| 417 | // alpha COOH and one for a residual chain chemical group, like | ||
| 418 | // epsilon NH2 for lysine. | ||
| 419 | |||
| 420 | ✗ | for(int kter = 0; kter < (*monomer_iterator_cst)->propList().size(); | |
| 421 | ++kter) | ||
| 422 | { | ||
| 423 | ✗ | Prop *prop = (*monomer_iterator_cst)->propList().at(kter); | |
| 424 | |||
| 425 | ✗ | if(prop->name() != "CHEMICAL_GROUP") | |
| 426 | ✗ | continue; | |
| 427 | |||
| 428 | // qDebug() << __FILE__ << __LINE__ | ||
| 429 | // << "Monomer has property CHEMICAL_GROUP..."; | ||
| 430 | |||
| 431 | // Get the chemical group out of the property. | ||
| 432 | |||
| 433 | ✗ | ChemicalGroup *chemicalGroup = | |
| 434 | ✗ | static_cast<ChemicalGroup *>(prop->data()); | |
| 435 | |||
| 436 | ✗ | if(static_cast<int>(chemicalGroup->getPolRule()) & | |
| 437 | static_cast<int>(Enums::ChemicalGroupTrapped::LEFT)) | ||
| 438 | { | ||
| 439 | // qDebug() << __FILE__ << __LINE__ | ||
| 440 | // << "... that is CHEMGROUP_LEFT_TRAPPED"; | ||
| 441 | |||
| 442 | // The chemical group we are dealing with is trapped | ||
| 443 | // when the monomer is polymerized on the left end, that | ||
| 444 | // is if the monomer is not the left end monomer of the | ||
| 445 | // sequence being analyzed. | ||
| 446 | |||
| 447 | // Thus we only can take it into account if one of | ||
| 448 | // two conditions are met: | ||
| 449 | |||
| 450 | // 1. The monomer is the left end monomer of the | ||
| 451 | // whole polymer sequence. | ||
| 452 | |||
| 453 | // 2. The monomer is the left end monomer of the | ||
| 454 | // region selection AND the selection type is | ||
| 455 | // oligomers(thus it does not get polymerized to | ||
| 456 | // the previous selection region). | ||
| 457 | |||
| 458 | ✗ | if(jter > 0) | |
| 459 | { | ||
| 460 | // Clearly we are not dealing with the left | ||
| 461 | // end of the polymer, so check if we have to | ||
| 462 | // account for this chemical group or not. | ||
| 463 | |||
| 464 | ✗ | if(!is_left_most_sequence_range) | |
| 465 | { | ||
| 466 | // The current libXpertMassCore::Coordinates is not the | ||
| 467 | // left-most libXpertMassCore::Coordinates in the | ||
| 468 | // libXpertMassCore::CoordinateList, thus we cannot | ||
| 469 | // consider it to be the "left end | ||
| 470 | // iter_index_range" of the | ||
| 471 | // libXpertMassCore::CoordinateList. Just continue without | ||
| 472 | // exploring any more. | ||
| 473 | ✗ | continue; | |
| 474 | } | ||
| 475 | ✗ | if(jter == start_index) | |
| 476 | { | ||
| 477 | // The current monomer is the first | ||
| 478 | // monomer of libXpertMassCore::Coordinates. We only take | ||
| 479 | // into account the chemical group if each | ||
| 480 | // libXpertMassCore::Coordinates is to be considered an | ||
| 481 | // oligomer. | ||
| 482 | |||
| 483 | ✗ | if(m_calcOptions.getSelectionType() != | |
| 484 | Enums::SelectionType::OLIGOMERS) | ||
| 485 | ✗ | continue; | |
| 486 | } | ||
| 487 | } | ||
| 488 | } | ||
| 489 | |||
| 490 | ✗ | if(static_cast<int>(chemicalGroup->getPolRule()) & | |
| 491 | static_cast<int>(Enums::ChemicalGroupTrapped::RIGHT)) | ||
| 492 | { | ||
| 493 | // qDebug() << __FILE__ << __LINE__ | ||
| 494 | // << "... that is CHEMGROUP_RIGHT_TRAPPED"; | ||
| 495 | |||
| 496 | // See explanations above. | ||
| 497 | |||
| 498 | ✗ | if(jter < (qsizetype) polymerSize - 1) | |
| 499 | { | ||
| 500 | // Clearly, we are not dealing with the right | ||
| 501 | // end of the polymer. | ||
| 502 | |||
| 503 | ✗ | if(!is_right_most_sequence_range) | |
| 504 | { | ||
| 505 | // The current libXpertMassCore::Coordinates is not the | ||
| 506 | // right-most libXpertMassCore::Coordinates of the | ||
| 507 | // libXpertMassCore::CoordinateList, thus we cannot | ||
| 508 | // consider it to be the "right end | ||
| 509 | // iter_index_range" of the | ||
| 510 | // libXpertMassCore::CoordinateList. Just continue without | ||
| 511 | // exploring anymore. | ||
| 512 | ✗ | continue; | |
| 513 | } | ||
| 514 | ✗ | if(jter == stop_index) | |
| 515 | { | ||
| 516 | // The current monomer is the last monomer | ||
| 517 | // of libXpertMassCore::Coordinates. We only take into | ||
| 518 | // account the chemical group if each | ||
| 519 | // libXpertMassCore::Coordinates is to be considered an | ||
| 520 | // oligomer(and not a residual chain). | ||
| 521 | |||
| 522 | ✗ | if(m_calcOptions.getSelectionType() != | |
| 523 | Enums::SelectionType::OLIGOMERS) | ||
| 524 | ✗ | continue; | |
| 525 | } | ||
| 526 | } | ||
| 527 | } | ||
| 528 | |||
| 529 | ✗ | if(iter == 0 && | |
| 530 | ✗ | static_cast<int>(m_calcOptions.getPolymerEntities()) & | |
| 531 | static_cast<int>(Enums::ChemicalEntity::LEFT_END_MODIF)) | ||
| 532 | { | ||
| 533 | // We are iterating in the monomer that is at the | ||
| 534 | // beginning of the polymer sequence. We should | ||
| 535 | // test if the chemical group we are dealing with | ||
| 536 | // right now has a special rule for the left end | ||
| 537 | // of the polymer sequence region. | ||
| 538 | |||
| 539 | ✗ | int ret = accountPolymerEndModif( | |
| 540 | Enums::ChemicalEntity::LEFT_END_MODIF, *chemicalGroup); | ||
| 541 | ✗ | if(ret >= 0) | |
| 542 | { | ||
| 543 | // qDebug() << __FILE__ << __LINE__ | ||
| 544 | // << "Accounted for left end modif."; | ||
| 545 | |||
| 546 | ✗ | processedChemicalGroups += ret; | |
| 547 | ✗ | continue; | |
| 548 | } | ||
| 549 | } | ||
| 550 | |||
| 551 | ✗ | if(iter == (qsizetype) polymerSize - 1 && | |
| 552 | ✗ | static_cast<int>(m_calcOptions.getPolymerEntities()) & | |
| 553 | static_cast<int>(Enums::ChemicalEntity::RIGHT_END_MODIF)) | ||
| 554 | { | ||
| 555 | ✗ | int ret = accountPolymerEndModif( | |
| 556 | Enums::ChemicalEntity::RIGHT_END_MODIF, *chemicalGroup); | ||
| 557 | ✗ | if(ret >= 0) | |
| 558 | { | ||
| 559 | // qDebug() << __FILE__ << __LINE__ | ||
| 560 | // << "Accounted for right end modif."; | ||
| 561 | |||
| 562 | ✗ | processedChemicalGroups += ret; | |
| 563 | ✗ | continue; | |
| 564 | } | ||
| 565 | } | ||
| 566 | |||
| 567 | ✗ | if(static_cast<int>(m_calcOptions.getMonomerEntities()) & | |
| 568 | ✗ | static_cast<int>(Enums::ChemicalEntity::MODIF) && | |
| 569 | ✗ | (*monomer_iterator_cst)->isModified()) | |
| 570 | { | ||
| 571 | ✗ | int ret = accountMonomerModif(*(*monomer_iterator_cst), | |
| 572 | *chemicalGroup); | ||
| 573 | ✗ | if(ret >= 0) | |
| 574 | { | ||
| 575 | // qDebug() << __FILE__ << __LINE__ | ||
| 576 | // << "Accounted for monomer modif."; | ||
| 577 | |||
| 578 | ✗ | processedChemicalGroups += ret; | |
| 579 | ✗ | continue; | |
| 580 | } | ||
| 581 | } | ||
| 582 | |||
| 583 | ✗ | double charge = calculateChargeRatio( | |
| 584 | ✗ | chemicalGroup->getPka(), chemicalGroup->isAcidCharged()); | |
| 585 | |||
| 586 | // qDebug() << __FILE__ << __LINE__ | ||
| 587 | // << "Charge:" << charge; | ||
| 588 | |||
| 589 | ✗ | if(charge < 0) | |
| 590 | ✗ | m_negativeCharges += charge; | |
| 591 | ✗ | else if(charge > 0) | |
| 592 | ✗ | m_positiveCharges += charge; | |
| 593 | |||
| 594 | // qDebug() << __FILE__ << __LINE__ | ||
| 595 | // << "Pos =" << m_positiveCharges | ||
| 596 | // << "Neg = " << m_negativeCharges; | ||
| 597 | |||
| 598 | ✗ | ++processedChemicalGroups; | |
| 599 | } | ||
| 600 | // End of | ||
| 601 | // for (int kter = 0; kter < monomer->propList().size(); ++kter) | ||
| 602 | |||
| 603 | // qDebug() << __FILE__ << __LINE__ | ||
| 604 | // << "End dealing with Monomer:" << | ||
| 605 | // seqMonomer->name() | ||
| 606 | // << "position:" << jter + 1; | ||
| 607 | ✗ | } | |
| 608 | // End of | ||
| 609 | // for (int jter = start_index; jter < stop_index + 1; ++jter) | ||
| 610 | |||
| 611 | // qDebug() << __FILE__ << __LINE__ | ||
| 612 | // << "End dealing with libXpertMassCore::Coordinates"; | ||
| 613 | } | ||
| 614 | // End of | ||
| 615 | // for (int iter = 0; iter < index_range_collection.size(); ++iter) | ||
| 616 | |||
| 617 | // We have finished processing all the libXpertMassCore::Coordinates in the list. | ||
| 618 | |||
| 619 | return processedChemicalGroups; | ||
| 620 | } | ||
| 621 | |||
| 622 | /*! | ||
| 623 | \brief Calculates the isoelectric point. | ||
| 624 | |||
| 625 | The isoelectric point is the pH at which a given analyte will have a net charge | ||
| 626 | of 0, that is, when the count of negative charges will be equal to the count of | ||
| 627 | positive charges. | ||
| 628 | |||
| 629 | The pI will be stored in the \l m_pi member variable. | ||
| 630 | |||
| 631 | Returns the count of chemical groups that have been processed. | ||
| 632 | |||
| 633 | \sa calculateCharges() | ||
| 634 | */ | ||
| 635 | int | ||
| 636 | ✗ | PkaPhPi::calculatePi() | |
| 637 | { | ||
| 638 | ✗ | int processedChemicalGroups = 0; | |
| 639 | ✗ | int iteration = 0; | |
| 640 | |||
| 641 | ✗ | double netCharge = 0; | |
| 642 | ✗ | double firstCharge = 0; | |
| 643 | ✗ | double thirdCharge = 0; | |
| 644 | |||
| 645 | // We of course need monomers ! Instead, we may not need modifs. | ||
| 646 | ✗ | if(!m_monomers.size()) | |
| 647 | { | ||
| 648 | ✗ | m_pi = 0; | |
| 649 | ✗ | m_ph = 0; | |
| 650 | |||
| 651 | ✗ | return -1; | |
| 652 | } | ||
| 653 | |||
| 654 | ✗ | m_ph = 0; | |
| 655 | |||
| 656 | ✗ | while(true) | |
| 657 | { | ||
| 658 | // qDebug() << "Current pH being tested:" << m_ph; | ||
| 659 | |||
| 660 | ✗ | processedChemicalGroups = calculateCharges(); | |
| 661 | |||
| 662 | ✗ | if(processedChemicalGroups == -1) | |
| 663 | { | ||
| 664 | ✗ | qDebug() << "Failed to calculate net charge for pH" << m_ph; | |
| 665 | |||
| 666 | ✗ | m_pi = 0; | |
| 667 | ✗ | m_ph = 0; | |
| 668 | |||
| 669 | ✗ | return -1; | |
| 670 | } | ||
| 671 | |||
| 672 | ✗ | netCharge = m_positiveCharges + m_negativeCharges; | |
| 673 | |||
| 674 | // Note that if the 0.01 tested_ph step is enough to switch the | ||
| 675 | // net charge from one excess value to another excess value in | ||
| 676 | // the opposite direction, we'll enter an infinite loop. | ||
| 677 | // | ||
| 678 | // The evidence for such loop is that every other two measures, | ||
| 679 | // the net charge of the polymer sequence will be the same. | ||
| 680 | // | ||
| 681 | // Here we test this so that we can break the loop. | ||
| 682 | |||
| 683 | |||
| 684 | ✗ | ++iteration; | |
| 685 | |||
| 686 | ✗ | if(iteration == 1) | |
| 687 | { | ||
| 688 | firstCharge = netCharge; | ||
| 689 | } | ||
| 690 | ✗ | else if(iteration == 3) | |
| 691 | { | ||
| 692 | ✗ | thirdCharge = netCharge; | |
| 693 | |||
| 694 | ✗ | if(firstCharge == thirdCharge) | |
| 695 | break; | ||
| 696 | |||
| 697 | iteration = 0; | ||
| 698 | |||
| 699 | firstCharge = netCharge; | ||
| 700 | } | ||
| 701 | |||
| 702 | // At this point we have to test the net charge: | ||
| 703 | |||
| 704 | ✗ | if(netCharge >= -0.1 && netCharge <= 0.1) | |
| 705 | { | ||
| 706 | // qDebug() << "Breaking loop with netCharge:" << netCharge; | ||
| 707 | |||
| 708 | break; | ||
| 709 | } | ||
| 710 | |||
| 711 | ✗ | if(netCharge > 0) | |
| 712 | { | ||
| 713 | // The test ph is too low. | ||
| 714 | |||
| 715 | ✗ | m_ph += 0.01; | |
| 716 | // qDebug() << "Set new pH m_ph += 0.01:" << m_ph | ||
| 717 | // << "netCharge:" << netCharge; | ||
| 718 | |||
| 719 | ✗ | continue; | |
| 720 | } | ||
| 721 | |||
| 722 | ✗ | if(netCharge < 0) | |
| 723 | { | ||
| 724 | // The test ph is too high. | ||
| 725 | |||
| 726 | ✗ | m_ph -= 0.01; | |
| 727 | // qDebug() << "Set new pH m_ph -= 0.01:" << m_ph | ||
| 728 | // << "netCharge:" << netCharge; | ||
| 729 | |||
| 730 | ✗ | continue; | |
| 731 | } | ||
| 732 | } | ||
| 733 | // End of | ||
| 734 | // while(true) | ||
| 735 | |||
| 736 | // At this point m_pi is m_ph. | ||
| 737 | |||
| 738 | ✗ | m_pi = m_ph; | |
| 739 | // qDebug() << "pi is:" << m_pi; | ||
| 740 | |||
| 741 | |||
| 742 | ✗ | return processedChemicalGroups; | |
| 743 | } | ||
| 744 | |||
| 745 | /*! | ||
| 746 | \brief Returns the ratio between the charged and the uncharged forms of the | ||
| 747 | chemical entity using \a pka and \a is_acid_charged. If the charge is negative, | ||
| 748 | the returned ratio is negative, positive otherwise. | ||
| 749 | |||
| 750 | The charged and uncharged species are the AH an A- species of the acido-basic | ||
| 751 | theory. | ||
| 752 | |||
| 753 | The calculation is based on the use of the m_ph member variable value. | ||
| 754 | */ | ||
| 755 | double | ||
| 756 | ✗ | PkaPhPi::calculateChargeRatio(double pka, bool is_acid_charged) | |
| 757 | { | ||
| 758 | ✗ | double aOverAh = 0; | |
| 759 | |||
| 760 | ✗ | if(pka < 0) | |
| 761 | return 0; | ||
| 762 | ✗ | if(pka > 14) | |
| 763 | return 0; | ||
| 764 | |||
| 765 | ✗ | if(m_ph < 0) | |
| 766 | return 0; | ||
| 767 | ✗ | if(m_ph > 14) | |
| 768 | return 0; | ||
| 769 | |||
| 770 | |||
| 771 | // Example with pKa = 4.25(Glu) ; pH = 4.16, thus we are more | ||
| 772 | // acidic than pKa, we expect AH to be greater than A by a small | ||
| 773 | // margin. | ||
| 774 | |||
| 775 | ✗ | aOverAh = (double)pow(10, (m_ph - pka)); | |
| 776 | // aOverAh = 0.81283051616409951 (confirmed manually) | ||
| 777 | |||
| 778 | ✗ | if(aOverAh < 1) | |
| 779 | { | ||
| 780 | /* The solution contains more acid forms(AH) than basic forms | ||
| 781 | (A). | ||
| 782 | */ | ||
| 783 | ✗ | if(is_acid_charged) | |
| 784 | ✗ | return (1 - aOverAh); | |
| 785 | else | ||
| 786 | // The acid is not charged, that is, it is a COOH. | ||
| 787 | // AH = 1 - A | ||
| 788 | // A = aOverAh.AH | ||
| 789 | // A = aOverAh.(1-A) | ||
| 790 | // A = aOverAh - aOverAh.A | ||
| 791 | // A(1+aOverAh) = aOverAh | ||
| 792 | // A = aOverAh /(1+aOverAh), for us this is | ||
| 793 | // A = 0.81283 / 1.81283 = 0.448 | ||
| 794 | |||
| 795 | // And not - aOverAh, that is - aOverAh. | ||
| 796 | |||
| 797 | // Below seems faulty(20071204) | ||
| 798 | // return(- aOverAh); | ||
| 799 | |||
| 800 | // Tentative correction(20071204) | ||
| 801 | ✗ | return (-(aOverAh / (1 + aOverAh))); | |
| 802 | } | ||
| 803 | ✗ | else if(aOverAh > 1) | |
| 804 | { | ||
| 805 | /* The solution contains more basic forms(A) than acid forms | ||
| 806 | (AH). | ||
| 807 | */ | ||
| 808 | ✗ | if(is_acid_charged) | |
| 809 | ✗ | return (1 / aOverAh); | |
| 810 | else | ||
| 811 | ✗ | return (-(1 - (1 / aOverAh))); | |
| 812 | } | ||
| 813 | ✗ | else if(aOverAh == 1) | |
| 814 | { | ||
| 815 | /* The solution contains as many acid forms(AH) as basic forms | ||
| 816 | (H). | ||
| 817 | */ | ||
| 818 | ✗ | if(is_acid_charged) | |
| 819 | return (aOverAh / 2); | ||
| 820 | else | ||
| 821 | ✗ | return (-aOverAh / 2); | |
| 822 | } | ||
| 823 | else | ||
| 824 | ✗ | qFatal("Programming error."); | |
| 825 | |||
| 826 | return 0; | ||
| 827 | } | ||
| 828 | |||
| 829 | /*! | ||
| 830 | \brief Accounts for the \a chemical_group according the the polymer chemical | ||
| 831 | entities defined in \a polymer_chem_ent. | ||
| 832 | |||
| 833 | A chemical group is described as follows: | ||
| 834 | |||
| 835 | \code | ||
| 836 | <monomer> | ||
| 837 | <code>C</code> | ||
| 838 | <mnmchemgroup> | ||
| 839 | <name>N-term NH2</name> | ||
| 840 | <pka>9.6</pka> | ||
| 841 | <acidcharged>TRUE</acidcharged> | ||
| 842 | <polrule>left_trapped</polrule> | ||
| 843 | <chemgrouprule> | ||
| 844 | <entity>LE_PLM_MODIF</entity> | ||
| 845 | <name>Acetylation</name> | ||
| 846 | <outcome>LOST</outcome> | ||
| 847 | </chemgrouprule> | ||
| 848 | </mnmchemgroup> | ||
| 849 | <mnmchemgroup> | ||
| 850 | <name>C-term COOH</name> | ||
| 851 | <pka>2.35</pka> | ||
| 852 | <acidcharged>FALSE</acidcharged> | ||
| 853 | <polrule>right_trapped</polrule> | ||
| 854 | </mnmchemgroup> | ||
| 855 | <mnmchemgroup> | ||
| 856 | <name>Lateral SH2</name> | ||
| 857 | <pka>8.3</pka> | ||
| 858 | <acidcharged>FALSE</acidcharged> | ||
| 859 | <polrule>never_trapped</polrule> | ||
| 860 | </mnmchemgroup> | ||
| 861 | </monomer> | ||
| 862 | \endcode | ||
| 863 | |||
| 864 | Returns the count of rules that were accounted for, or -1 if none was. | ||
| 865 | */ | ||
| 866 | int | ||
| 867 | ✗ | PkaPhPi::accountPolymerEndModif(Enums::ChemicalEntity polymer_chem_ent, | |
| 868 | const ChemicalGroup &chemical_group) | ||
| 869 | { | ||
| 870 | ✗ | QString modifName; | |
| 871 | ✗ | ChemicalGroupRuleSPtr chemical_group_rule_sp = nullptr; | |
| 872 | ✗ | int count = 0; | |
| 873 | |||
| 874 | // Get the name of the modification of the polymer (if any) and get | ||
| 875 | // the rule dealing with that polymer modification (if any). | ||
| 876 | |||
| 877 | ✗ | std::size_t chem_group_rule_index; | |
| 878 | |||
| 879 | ✗ | if(static_cast<int>(polymer_chem_ent) & | |
| 880 | static_cast<int>(Enums::ChemicalEntity::LEFT_END_MODIF)) | ||
| 881 | { | ||
| 882 | ✗ | modifName = mcsp_polymer->getLeftEndModifCstRef().getName(); | |
| 883 | ✗ | chemical_group_rule_sp = chemical_group.findRuleByEntityAndName( | |
| 884 | ✗ | "LE_PLM_MODIF", modifName, chem_group_rule_index); | |
| 885 | } | ||
| 886 | ✗ | if(static_cast<int>(polymer_chem_ent) & | |
| 887 | static_cast<int>(Enums::ChemicalEntity::RIGHT_END_MODIF)) | ||
| 888 | { | ||
| 889 | ✗ | modifName = mcsp_polymer->getRightEndModifCstRef().getName(); | |
| 890 | ✗ | chemical_group_rule_sp = chemical_group.findRuleByEntityAndName( | |
| 891 | ✗ | "RE_PLM_MODIF", modifName, chem_group_rule_index); | |
| 892 | } | ||
| 893 | else | ||
| 894 | ✗ | qFatal("Programming error."); | |
| 895 | |||
| 896 | |||
| 897 | // The polymer might not be modified, and also the chemical group | ||
| 898 | // passed as parameter might not contain any rule about any polymer | ||
| 899 | // modification. In that case we just have nothing to do. | ||
| 900 | |||
| 901 | ✗ | if(modifName.isEmpty()) | |
| 902 | { | ||
| 903 | ✗ | if(chemical_group_rule_sp != nullptr) | |
| 904 | { | ||
| 905 | ✗ | double charge = calculateChargeRatio(chemical_group.getPka(), | |
| 906 | ✗ | chemical_group.isAcidCharged()); | |
| 907 | ✗ | if(charge < 0) | |
| 908 | ✗ | m_negativeCharges += charge; | |
| 909 | ✗ | else if(charge > 0) | |
| 910 | ✗ | m_positiveCharges += charge; | |
| 911 | |||
| 912 | return ++count; | ||
| 913 | } | ||
| 914 | else | ||
| 915 | { | ||
| 916 | // The polymer end was NOT modified and the chemical group | ||
| 917 | // was NOT eligible for a polymer end modification. This | ||
| 918 | // means that we do not have to process it, and we return -1 | ||
| 919 | // so that the caller function knows we did not do anything | ||
| 920 | // and that this chemical group should continue to undergo | ||
| 921 | // analysis without skipping it. | ||
| 922 | |||
| 923 | return -1; | ||
| 924 | } | ||
| 925 | } | ||
| 926 | // End of | ||
| 927 | // if (modifName.isEmpty()) | ||
| 928 | |||
| 929 | ✗ | if(chemical_group_rule_sp == nullptr) | |
| 930 | { | ||
| 931 | // This chemical group was not "designed" to receive any polymer | ||
| 932 | // end modification, so we have nothing to do with it and we | ||
| 933 | // return -1 so that the caller function knows we did not do | ||
| 934 | // anything and that this chemical group should continue to | ||
| 935 | // undergo analysis without skipping it. | ||
| 936 | |||
| 937 | return -1; | ||
| 938 | } | ||
| 939 | |||
| 940 | // At this point we know that the chemical group 'group' we are | ||
| 941 | // analyzing is eligible for a polymer left end modification and | ||
| 942 | // that it is indeed modified with a correcct modification. So we | ||
| 943 | // have a rule for it. Let's continue the analysis. | ||
| 944 | |||
| 945 | // Apparently the rule has data matching the ones we are looking | ||
| 946 | // for. At this point we should now what action to take for this | ||
| 947 | // group. | ||
| 948 | |||
| 949 | ✗ | if(chemical_group_rule_sp->getFate() == Enums::ChemicalGroupFate::LOST) | |
| 950 | { | ||
| 951 | // We do not use the current chemical group 'group' because the | ||
| 952 | // polymer end's modification has abolished it. | ||
| 953 | } | ||
| 954 | ✗ | else if(chemical_group_rule_sp->getFate() == Enums::ChemicalGroupFate::PRESERVED) | |
| 955 | { | ||
| 956 | ✗ | double charge = calculateChargeRatio(chemical_group.getPka(), | |
| 957 | ✗ | chemical_group.isAcidCharged()); | |
| 958 | ✗ | if(charge < 0) | |
| 959 | ✗ | m_negativeCharges += charge; | |
| 960 | ✗ | else if(charge > 0) | |
| 961 | ✗ | m_positiveCharges += charge; | |
| 962 | |||
| 963 | ✗ | return ++count; | |
| 964 | } | ||
| 965 | else | ||
| 966 | ✗ | qFatal("Programming error."); | |
| 967 | |||
| 968 | // Whatever we should do with the left/right end monomer's chemgroup, | ||
| 969 | // we should take into account the modification itself that might | ||
| 970 | // have brought chemgroup(s) worth calculating their intrinsic | ||
| 971 | // charges! | ||
| 972 | |||
| 973 | // Find a modif object in the local list of modif objects, that has | ||
| 974 | // the same name as the modification with which the left/right end | ||
| 975 | // of the polymer is modified. We'll see what chemgroup(s) this | ||
| 976 | // modification brings to the polymer sequence. | ||
| 977 | |||
| 978 | ✗ | std::vector<ModifCstSPtr>::const_iterator modif_iterator_cst = | |
| 979 | ✗ | std::find_if(m_modifs.cbegin(), | |
| 980 | m_modifs.cend(), | ||
| 981 | ✗ | [modifName](const ModifCstSPtr &modif_csp) { | |
| 982 | ✗ | return modif_csp->getName() == modifName; | |
| 983 | }); | ||
| 984 | |||
| 985 | ✗ | if(modif_iterator_cst == m_modifs.cend()) | |
| 986 | { | ||
| 987 | // qDebug() << __FILE__ << __LINE__ | ||
| 988 | // << "Information: following modif not in local list:" | ||
| 989 | // << modifName; | ||
| 990 | |||
| 991 | return count; | ||
| 992 | } | ||
| 993 | |||
| 994 | ✗ | for(int jter = 0; jter < (*modif_iterator_cst)->propList().size(); ++jter) | |
| 995 | { | ||
| 996 | ✗ | Prop *prop = (*modif_iterator_cst)->propList().at(jter); | |
| 997 | |||
| 998 | ✗ | if(prop->name() != "CHEMICAL_GROUP") | |
| 999 | ✗ | continue; | |
| 1000 | |||
| 1001 | // Get the chemical group out of the property. | ||
| 1002 | |||
| 1003 | ✗ | const ChemicalGroup *chemicalGroup = | |
| 1004 | ✗ | static_cast<const ChemicalGroup *>(prop->data()); | |
| 1005 | |||
| 1006 | ✗ | double charge = calculateChargeRatio(chemicalGroup->getPka(), | |
| 1007 | ✗ | chemicalGroup->isAcidCharged()); | |
| 1008 | ✗ | if(charge < 0) | |
| 1009 | ✗ | m_negativeCharges += charge; | |
| 1010 | ✗ | else if(charge > 0) | |
| 1011 | ✗ | m_positiveCharges += charge; | |
| 1012 | |||
| 1013 | ✗ | ++count; | |
| 1014 | } | ||
| 1015 | |||
| 1016 | return count; | ||
| 1017 | ✗ | } | |
| 1018 | |||
| 1019 | |||
| 1020 | /*! | ||
| 1021 | \brief Accounts for the \a chemical_group in the context of the \a monomer. | ||
| 1022 | |||
| 1023 | A chemical group is described as follows: | ||
| 1024 | |||
| 1025 | \code | ||
| 1026 | <monomer> | ||
| 1027 | <code>C</code> | ||
| 1028 | <mnmchemgroup> | ||
| 1029 | <name>N-term NH2</name> | ||
| 1030 | <pka>9.6</pka> | ||
| 1031 | <acidcharged>TRUE</acidcharged> | ||
| 1032 | <polrule>left_trapped</polrule> | ||
| 1033 | <chemgrouprule> | ||
| 1034 | <entity>LE_PLM_MODIF</entity> | ||
| 1035 | <name>Acetylation</name> | ||
| 1036 | <outcome>LOST</outcome> | ||
| 1037 | </chemgrouprule> | ||
| 1038 | </mnmchemgroup> | ||
| 1039 | <mnmchemgroup> | ||
| 1040 | <name>C-term COOH</name> | ||
| 1041 | <pka>2.35</pka> | ||
| 1042 | <acidcharged>FALSE</acidcharged> | ||
| 1043 | <polrule>right_trapped</polrule> | ||
| 1044 | </mnmchemgroup> | ||
| 1045 | <mnmchemgroup> | ||
| 1046 | <name>Lateral SH2</name> | ||
| 1047 | <pka>8.3</pka> | ||
| 1048 | <acidcharged>FALSE</acidcharged> | ||
| 1049 | <polrule>never_trapped</polrule> | ||
| 1050 | </mnmchemgroup> | ||
| 1051 | </monomer> | ||
| 1052 | \endcode | ||
| 1053 | |||
| 1054 | Returns the count of rules that were accounted for, or -1 if none was. | ||
| 1055 | */ | ||
| 1056 | int | ||
| 1057 | ✗ | PkaPhPi::accountMonomerModif(const Monomer &monomer, | |
| 1058 | ChemicalGroup &chemical_group) | ||
| 1059 | { | ||
| 1060 | ✗ | QString modifName; | |
| 1061 | ✗ | ChemicalGroupRuleSPtr rule_sp = nullptr; | |
| 1062 | |||
| 1063 | ✗ | int count = 0; | |
| 1064 | |||
| 1065 | // For each modification in the monomer, make the accounting work. | ||
| 1066 | |||
| 1067 | ✗ | std::size_t chem_group_rule_index = 0; | |
| 1068 | |||
| 1069 | ✗ | for(const ModifSPtr &modif_sp : monomer.getModifsCstRef()) | |
| 1070 | { | ||
| 1071 | // Get the name of the modification of the monomer(if any) and get | ||
| 1072 | // the rule dealing with that monomer modification(if any). | ||
| 1073 | |||
| 1074 | ✗ | modifName = modif_sp->getName(); | |
| 1075 | |||
| 1076 | ✗ | rule_sp = chemical_group.findRuleByEntityAndName( | |
| 1077 | ✗ | "MONOMER_MODIF", modifName, chem_group_rule_index); | |
| 1078 | |||
| 1079 | ✗ | if(modifName.isEmpty()) | |
| 1080 | { | ||
| 1081 | // The monomer does not seem to be modified. However, we still | ||
| 1082 | // have to make sure that the chemgroup that we were parsing is | ||
| 1083 | // actually a chemgroup suitable for a modif. If this chemgroup | ||
| 1084 | // was actually suitable for a monomer modif, but it is not | ||
| 1085 | // effectively modified, that means that we have to calculate | ||
| 1086 | // the charge for the non-modified chemgroup... | ||
| 1087 | |||
| 1088 | ✗ | if(rule_sp != nullptr) | |
| 1089 | { | ||
| 1090 | ✗ | double charge = calculateChargeRatio( | |
| 1091 | ✗ | chemical_group.getPka(), chemical_group.isAcidCharged()); | |
| 1092 | ✗ | if(charge < 0) | |
| 1093 | ✗ | m_negativeCharges += charge; | |
| 1094 | ✗ | else if(charge > 0) | |
| 1095 | ✗ | m_positiveCharges += charge; | |
| 1096 | |||
| 1097 | ✗ | return ++count; | |
| 1098 | } | ||
| 1099 | else | ||
| 1100 | { | ||
| 1101 | // The current monomer was NOT modified, and the chemgroup | ||
| 1102 | // was NOT eligible for a monomer modification. This means | ||
| 1103 | // that we do not have to process it, and we return -1 so | ||
| 1104 | // that the caller function knows we did not do anything and | ||
| 1105 | // that this chemgroup should continue to undergo analysis | ||
| 1106 | // without skipping it. | ||
| 1107 | |||
| 1108 | return -1; | ||
| 1109 | } | ||
| 1110 | } | ||
| 1111 | // End of | ||
| 1112 | // if (modifName.isEmpty()) | ||
| 1113 | |||
| 1114 | ✗ | if(rule_sp == nullptr) | |
| 1115 | { | ||
| 1116 | // This chemgroup was not "designed" to receive any | ||
| 1117 | // modification, so we have nothing to do with it, and we return | ||
| 1118 | // -1 to let the caller know that its processing should be | ||
| 1119 | // continued in the caller's function space. | ||
| 1120 | |||
| 1121 | return -1; | ||
| 1122 | } | ||
| 1123 | |||
| 1124 | // At this point, we know that the chemgroup we are analyzing is | ||
| 1125 | // eligible for a modification and that we have a rule for it. Let's | ||
| 1126 | // continue the analysis: | ||
| 1127 | |||
| 1128 | // Apparently, a rule object has member data matching the ones we | ||
| 1129 | // were looking for. At this point we should know what action to | ||
| 1130 | // take for this chemgroup. | ||
| 1131 | |||
| 1132 | ✗ | if(rule_sp->getFate() == Enums::ChemicalGroupFate::LOST) | |
| 1133 | { | ||
| 1134 | // We do not use the current chemical group 'group' because the | ||
| 1135 | // monomer modification has abolished it. | ||
| 1136 | } | ||
| 1137 | ✗ | else if(rule_sp->getFate() == Enums::ChemicalGroupFate::PRESERVED) | |
| 1138 | { | ||
| 1139 | ✗ | double charge = calculateChargeRatio(chemical_group.getPka(), | |
| 1140 | ✗ | chemical_group.isAcidCharged()); | |
| 1141 | ✗ | if(charge < 0) | |
| 1142 | ✗ | m_negativeCharges += charge; | |
| 1143 | ✗ | else if(charge > 0) | |
| 1144 | ✗ | m_positiveCharges += charge; | |
| 1145 | |||
| 1146 | ✗ | return ++count; | |
| 1147 | } | ||
| 1148 | else | ||
| 1149 | ✗ | qFatalStream() << "Programming error."; | |
| 1150 | |||
| 1151 | // Whatever we should do with this monomer's chemgroup, we should | ||
| 1152 | // take into account the modification itself that might have brought | ||
| 1153 | // chemgroup(s) worth calculating their intrinsic charges! | ||
| 1154 | |||
| 1155 | // Find a modif object in the local list of modif objects, that has | ||
| 1156 | // the same name as the modification with which the monomer is | ||
| 1157 | // modified. We'll see what chemgroup(s) this modification brings to | ||
| 1158 | // the polymer sequence. | ||
| 1159 | |||
| 1160 | ✗ | std::vector<ModifCstSPtr>::const_iterator modif_monomer_iterator_cst = | |
| 1161 | ✗ | std::find_if(m_modifs.cbegin(), | |
| 1162 | m_modifs.cend(), | ||
| 1163 | ✗ | [modifName](const ModifCstSPtr &modif_csp) { | |
| 1164 | ✗ | return modif_csp->getName() == modifName; | |
| 1165 | }); | ||
| 1166 | |||
| 1167 | ✗ | if(modif_monomer_iterator_cst == m_modifs.cend()) | |
| 1168 | { | ||
| 1169 | // qDebug() << __FILE__ << __LINE__ | ||
| 1170 | // << "Information: following modif not in local list:" | ||
| 1171 | // << modifName; | ||
| 1172 | |||
| 1173 | return count; | ||
| 1174 | } | ||
| 1175 | |||
| 1176 | ✗ | for(int jter = 0; jter < (*modif_monomer_iterator_cst)->propList().size(); | |
| 1177 | ++jter) | ||
| 1178 | { | ||
| 1179 | ✗ | Prop *prop = (*modif_monomer_iterator_cst)->propList().at(jter); | |
| 1180 | |||
| 1181 | ✗ | if(prop->name() != "CHEMICAL_GROUP") | |
| 1182 | ✗ | continue; | |
| 1183 | |||
| 1184 | // Get the chemical group out of the property. | ||
| 1185 | |||
| 1186 | ✗ | const ChemicalGroup *chemicalGroup = | |
| 1187 | ✗ | static_cast<const ChemicalGroup *>(prop->data()); | |
| 1188 | |||
| 1189 | ✗ | double charge = calculateChargeRatio(chemicalGroup->getPka(), | |
| 1190 | ✗ | chemicalGroup->isAcidCharged()); | |
| 1191 | ✗ | if(charge < 0) | |
| 1192 | ✗ | m_negativeCharges += charge; | |
| 1193 | ✗ | else if(charge > 0) | |
| 1194 | ✗ | m_positiveCharges += charge; | |
| 1195 | |||
| 1196 | ✗ | ++count; | |
| 1197 | } | ||
| 1198 | } | ||
| 1199 | |||
| 1200 | return count; | ||
| 1201 | ✗ | } | |
| 1202 | |||
| 1203 | } // namespace libXpertMassCore | ||
| 1204 | } // namespace MsXpS | ||
| 1205 |