GCC Code Coverage Report


./
File: src/XpertMass/PkaPhPi.cpp
Date: 2024-08-24 11:26:06
Lines:
0/273
0.0%
Functions:
0/15
0.0%
Branches:
0/296
0.0%

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