massxpert-c229f4a1abde3c20b83a90e50f9c5d79104dfa5f/lib/oligomer.cpp

/* massXpert - the true massist's program.
   --------------------------------------
   Copyright(C) 2006,2007 Filippo Rusconi

   http://www.massxpert.org/massXpert

   This file is part of the massXpert project.

   The massxpert project is the successor to the "GNU polyxmass"
   project that is an official GNU project package(see
   www.gnu.org). The massXpert project is not endorsed by the GNU
   project, although it is released ---in its entirety--- under the
   GNU General Public License. A huge part of the code in massXpert
   is actually a C++ rewrite of code in GNU polyxmass. As such
   massXpert was started at the Centre National de la Recherche
   Scientifique(FRANCE), that granted me the formal authorization to
   publish it under this Free Software License.

   This software is free software; you can redistribute it and/or
   modify it under the terms of the GNU  General Public
   License version 3, as published by the Free Software Foundation.


   This software is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this software; if not, write to the

   Free Software Foundation, Inc.,

   51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
*/


/////////////////////// Local includes
#include "oligomer.hpp"
#include "polymer.hpp"
#include "ionizeRule.hpp"


namespace massXpert
{

  //! Constructs an oligomer.
  /*!

    \param polymer Polymer in which \c this oligomer spans a region. Cannot
    be 0.

    \param ponderable Ponderable used to initialized the data in \c
    this. Note that the initialization is shallow, as the data in the
    property list are not copied.

    \param name Name of the oligomer.

    \param startIndex Index of the first monomer of \c this oligomer in
    the polymer sequence(note that this index is stored in the very first
    Coordinates item in the CoordinateList).

    \param endIndex Index of the last monomer of \c this oligomer in the
    polymer sequence(note that this index is stored in the very first
    Coordinates item in the CoordinateList).

    \param ionizeRule IonizeRule to be used to initialize the Ionizable ancestor.

    \param isIonized Indicates if the oligomer should be considered
    ionized.
  */
  Oligomer::Oligomer(Polymer *polymer ,
                     const QString &name,
                     const QString &description,
                     const Ponderable &ponderable,
                     const IonizeRule &ionizeRule,
                     const CalcOptions & calcOptions,
                     bool isIonized,
                     int startIndex, int endIndex)
    : Ionizable(polymer->polChemDef(), name, ponderable,
                ionizeRule, isIonized),
      mp_polymer(polymer),
      m_description(description),
      m_calcOptions(calcOptions)
  {
    setStartEndIndices(startIndex, endIndex);
  }


  //! Constructs an oligomer.
  /*! The oligomer is constructed with its \c mp_polymer member set to
    0. This means that no reference to the polymer can be done. This
    constructor is useful when the oligomer must be created out of any
    polymer context, like when oligomers are created starting from raw
    text data in the MzLab window.

    \param polChemDef Polymer chemistry definition on which \c this
    oligomer is based. Cannot be 0.

    \param ponderable Ponderable used to initialized the data in \c
    this. Note that the initialization is shallow, as the data in the
    property list are not copied.

    \param name Name of the oligomer.

    \param startIndex Index of the first monomer of \c this oligomer in
    the polymer sequence(note that this index is stored in the very first
    Coordinates item in the CoordinateList).

    \param endIndex Index of the last monomer of \c this oligomer in the
    polymer sequence(note that this index is stored in the very first
    Coordinates item in the CoordinateList).

    \param ionizeRule IonizeRule to be used to initialize the Ionizable ancestor.

    \param isIonized Indicates if the oligomer should be considered
    ionized.
  */
  Oligomer::Oligomer(const PolChemDef *polChemDef,
                     const QString &name,
                     const QString &description,
                     const Ponderable &ponderable,
                     const IonizeRule &ionizeRule,
                     const CalcOptions &calcOptions,
                     bool isIonized,
                     int startIndex, int endIndex)
    : Ionizable(polChemDef, name, ponderable, ionizeRule, isIonized),
      mp_polymer(0),
      m_description(description),
      m_calcOptions(calcOptions)
  {
    setStartEndIndices(startIndex, endIndex);
    // if (startIndex < 0)
    //   qDebug() << __FILE__ << __LINE__
    //            << "Construct with startIndex:" << startIndex;
  }


  //! Constructs an oligomer.
  /*!

    \param polymer Polymer in which \c this oligomer spans a region. Cannot
    be 0.

    \param ponderable Ponderable used to initialized the data in \c
    this. Note that the initialization is shallow, as the data in the
    property list are not copied.

    \param name Name of the oligomer.

    \param startIndex Index of the first monomer of \c this oligomer in
    the polymer sequence(note that this index is stored in the very first
    Coordinates item in the CoordinateList).

    \param endIndex Index of the last monomer of \c this oligomer in the
    polymer sequence(note that this index is stored in the very first
    Coordinates item in the CoordinateList).
  */
  Oligomer::Oligomer(Polymer *polymer,
                     const QString &name,
                     const QString &description,
                     const Ponderable &ponderable,
                     int startIndex, int endIndex,
                     const CalcOptions &calcOptions)
    : Ionizable(polymer->polChemDef(), name, ponderable),
      mp_polymer(polymer),
      m_description(description),
      m_calcOptions(calcOptions)
  {
    setStartEndIndices(startIndex, endIndex);
    if (startIndex < 0)
      qDebug() << __FILE__ << __LINE__
               << "Construct with startIndex:" << startIndex;
  }


  //! Constructs an oligomer.
  /*! The oligomer is constructed with its \c mp_polymer member set to
    0. This means that no reference to the polymer can be done. This
    constructor is useful when the oligomer must be created out of any
    polymer context, like when oligomers are created starting from raw
    text data in the MzLab window.

    \param polChemDef Polymer chemistry definition on which \c this
    oligomer is based. Cannot be 0.

    \param ponderable Ponderable used to initialized the data in \c
    this. Note that the initialization is shallow, as the data in the
    property list are not copied.

    \param name Name of the oligomer.

    \param startIndex Index of the first monomer of \c this oligomer in
    the polymer sequence(note that this index is stored in the very first
    Coordinates item in the CoordinateList).

    \param endIndex Index of the last monomer of \c this oligomer in the
    polymer sequence(note that this index is stored in the very first
    Coordinates item in the CoordinateList).
  */
  Oligomer::Oligomer(const PolChemDef *polChemDef,
                     const QString &name,
                     const QString &description,
                     const Ponderable &ponderable,
                     const CalcOptions &calcOptions,
                     int startIndex, int endIndex)
    : Ionizable(polChemDef, name, ponderable),
      mp_polymer(0),
      m_description(description),
      m_calcOptions(calcOptions)
  {
    setStartEndIndices(startIndex, endIndex);
    if (startIndex < 0)
      qDebug() << __FILE__ << __LINE__
               << "Construct with startIndex:" << startIndex;
  }


  Oligomer::Oligomer(const Ionizable& ionizable,
                     const CalcOptions &calcOptions,
                     int startIndex, int endIndex)
    : Ionizable(ionizable),
      mp_polymer(0),
      m_description("NOT_SET"),
      m_calcOptions(calcOptions)
  {
    setStartEndIndices(startIndex, endIndex);
    if (startIndex < 0)
      qDebug() << __FILE__ << __LINE__
               << "Construct with startIndex:" << startIndex;
  }


  //! Constructs an oligomer.
  /*!

    \param polymer Polymer in which \c this oligomer spans a region. Cannot
    be 0.

    \param name Name of the oligomer.

    \param startIndex Index of the first monomer of \c this oligomer in
    the polymer sequence.

    \param endIndex Index of the last monomer of \c this oligomer in the
    polymer sequence.

    \param mono Monoisotopic mass.

    \param avg Average mass.
  */
  Oligomer::Oligomer(Polymer *polymer ,
                     const QString &name,
                     const QString &description,
                     double mono, double avg,
                     int startIndex, int endIndex,
                     const CalcOptions &calcOptions)
    : Ionizable(mp_polymer->polChemDef(), name, Ponderable(mono, avg)),
      mp_polymer(polymer),
      m_description(description),
      m_calcOptions(calcOptions)
  {
    setStartEndIndices(startIndex, endIndex);
    if (startIndex < 0)
      qDebug() << __FILE__ << __LINE__
               << "Construct with startIndex:" << startIndex;
  }


  //! Constructs an oligomer.
  /*! The oligomer is constructed with its \c mp_polymer member set to
    0. This means that no reference to the polymer can be done. This
    constructor is useful when the oligomer must be created out of any
    polymer context, like when oligomers are created starting from raw
    text data in the MzLab window.

    \param polChemDef Polymer chemistry definition on which \c this
    oligomer is based. Cannot be 0.

    \param name Name of the oligomer.

    \param startIndex Index of the first monomer of \c this oligomer in
    the polymer sequence.

    \param endIndex Index of the last monomer of \c this oligomer in the
    polymer sequence.

    \param mono Monoisotopic mass.

    \param avg Average mass.
  */
  Oligomer::Oligomer(const PolChemDef *polChemDef,
                     const QString &name,
                     const QString &description,
                     const CalcOptions &calcOptions,
                     double mono, double avg,
                     int startIndex, int endIndex)
    : Ionizable(polChemDef, name, Ponderable(mono, avg)),
      mp_polymer(0),
      m_description(description),
      m_calcOptions(calcOptions)
  {
    setStartEndIndices(startIndex, endIndex);
  }


  //! Constructs a copy of \p other.
  /*!  The copying is shallow, as the data in the property list are not
    copied.

    \param other oligomer to be used as a mold.
  */
  Oligomer::Oligomer(const Oligomer &other)
    : CoordinateList(other), Ionizable(other), PropListHolder(other),
      mp_polymer(other.mp_polymer),
      m_description(other.m_description),
      m_calcOptions(other.m_calcOptions)
  {
  }


  //! Destroys the oligomer.
  Oligomer::~Oligomer()
  {
    //   qDebug() << "~Oligomer()";
  }


  //! Returns the polymer.
  /*! \return the polymer.
   */
  const Polymer *
  Oligomer::polymer() const
  {
    return mp_polymer;
  }


  //! Sets the start and end indices.
  /*! \param value1 New start index
    \param value2 New end index
  */
  void
  Oligomer::setStartEndIndices(int value1, int value2)
  {
    if (!CoordinateList::size())
      {
	Coordinates *coordinates = new Coordinates(value1, value2);
	append(coordinates);

	//       qDebug() << __FILE__ << __LINE__
	// 		<< "[start--end]:" << startIndex() << endIndex();
      }
    else
      {
	Coordinates *coordinates = first();
	coordinates->setStart(value1);
	coordinates->setEnd(value2);

	//       qDebug() << __FILE__ << __LINE__
	// 		<< "[start--end]:" << startIndex() << endIndex();
      }
  }


  //! Sets the start index.
  /*! \param value New start index.
   */
  void
  Oligomer::setStartIndex(int value)
  {
    if (value < 0)
      qDebug() << __FILE__ << __LINE__
               << "setStartIndex with value:" << value;

    if (!CoordinateList::size())
      {
	Coordinates *coordinates = new Coordinates();
	coordinates->setStart(value);
	append(coordinates);
      }
    else
      {
	Coordinates *coordinates = first();
	coordinates->setStart(value);
      }
  }


  //! Returns the start index.
  /*! \return the start index.
   */
  int
  Oligomer::startIndex() const
  {
    if (!CoordinateList::size())
      return -1;

    Coordinates *coordinates = first();
    return coordinates->start();
  }


  //! Sets the end index.
  /*! \param value New end index.
   */
  void
  Oligomer::setEndIndex(int value)
  {
    if (!CoordinateList::size())
      {
	Coordinates *coordinates = new Coordinates();
	coordinates->setEnd(value);
	append(coordinates);
      }
    else
      {
	Coordinates *coordinates = first();
	coordinates->setEnd(value);
      }
  }


  //! Returns the end index.
  /*! \return the end index.
   */
  int
  Oligomer::endIndex() const
  {
    if (!CoordinateList::size())
      return -1;

    Coordinates *coordinates = first();
    return coordinates->end();
  }


  void
  Oligomer::setDescription(const QString &description)
  {
    m_description = description;
  }


  QString
  Oligomer::description() const
  {
    return m_description;
  }


  int
  Oligomer::appendCoordinates(CoordinateList *list)
  {
    Q_ASSERT(list);

    int count = 0;

    for (int iter = 0; iter < list->size(); ++iter)
      {
	Coordinates *iterCoordinates = list->at(iter);

	Coordinates *coordinates = new Coordinates(*iterCoordinates);

	append(coordinates);

	++count;
      }

    return count;
  }


  void
  Oligomer::setIonizeRule(IonizeRule &ionizeRule)
  {
    m_ionizeRule = ionizeRule;
  }


  IonizeRule &
  Oligomer::ionizeRule()
  {
    return m_ionizeRule;
  }

  void
  Oligomer::setCalcOptions(const CalcOptions &calcOptions)
  {
    m_calcOptions = calcOptions;
  }


  const CalcOptions &
  Oligomer::calcOptions() const
  {
    return m_calcOptions;
  }

  void
  Oligomer::updateCalcOptions()
  {
    // The data that need update are the CoordinateList elements
    // depending on the internal ::OligomerList data.

    m_calcOptions.setCoordinateList(*(static_cast<CoordinateList *>(this)));
  }


  const Monomer &
  Oligomer::atLeftEnd() const
  {
//     qDebug() << __FILE__ << __LINE__ << "Going to call at() with value"
// 	   << startIndex();

    return *(mp_polymer->at(startIndex()));
  }


  const Monomer &
  Oligomer::atRightEnd() const
  {
//     qDebug() << __FILE__ << __LINE__ << "Going to call at() with value"
// 	   << endIndex();

    return *(mp_polymer->at(endIndex()));


  }


  const Monomer *
  Oligomer::monomerAt(int index) const
  {
    Q_ASSERT(index >= 0);
    Q_ASSERT(index < mp_polymer->size() - 1);

//     qDebug() << __FILE__ << __LINE__ << "Going to call at() with value"
// 	      << index;

    return mp_polymer->at(index);
  }


  QList<CrossLink *> *
  Oligomer::crossLinkList()
  {
    return &m_crossLinkList;
  }


  bool
  Oligomer::addCrossLink(CrossLink *crossLink)
  {
    // Add the cross-link only if it does not exist already. Return true
    // only if the crossLink has been added.

    if (!m_crossLinkList.contains(crossLink))
      {
	m_crossLinkList.append(crossLink);

	return true;
      }

    return false;
  }


  // ELEMENTAL CALCULATION FUNCTION
  /////////////////////////////////

  QString
  Oligomer::elementalComposition()
  {
    int times = 0;

    if (m_calcOptions.selectionType() == SELECTION_TYPE_RESIDUAL_CHAINS)
      {
	times = 1;

        // qDebug() << __FILE__ << __LINE__
        //          << "SELECTION_TYPE_RESIDUAL_CHAINS ; times:" << times;
      }
    else
      {
        // 	times = CoordinateList::size();

        // Use the version whereby we can perform a sanity check that
        // m_calcOptions was updated with the proper CoordinateList
        // data.

        times = m_calcOptions.coordinateList().size();
        if(times != CoordinateList::size())
          qFatal("Fatal error at %s@%d. Aborting.",__FILE__, __LINE__);

        // qDebug() << __FILE__ << __LINE__
        //          << "SELECTION_TYPE_OLIGOMERS ; times:" << times;
      }

    return mp_polymer->elementalComposition
      (m_ionizeRule,
       *(static_cast<CoordinateList *>(this)),
       m_calcOptions);
  }


  /////////////////////////////////
  // ELEMENTAL CALCULATION FUNCTION


  int
  Oligomer::makeMonomerText()
  {
    // Prepare the text version of the oligomer's sequence by basing
    // it on the coordinates of *this oligomer and set that text to
    // m_monomerText.

    if (!mp_polymer)
      return 0;

    QString *text = monomerText();

    m_monomerText = *text;

    delete(text);

    return m_monomerText.size();
  }


  QString *
  Oligomer::monomerText()
  {
    // Allocate a new string to hold the text version of *this
    // oligomer's sequence.

    // There are two situations:

    // 1. The mp_polymer member is non-0, we can get access to it. Ask
    // the polymer to do the work. This is the most faithful
    // situation. But the caller must first ensure that the polymer
    // actually exists.

    // 2. The mp_polymer member is 0, we may have the oligomer
    // sequence stored in *this oligomer. Test that.

    QString *text = new QString();

    if (mp_polymer)
      {
	// For each oligomer(if more than one, which is the case when the
	// oligomer is actually a cross-linked oligomer), create a string...

	int oligomerCount = CoordinateList::size();

	for(int iter = 0; iter < oligomerCount; ++iter)
	  {
	    Coordinates *coordinates = CoordinateList::at(iter);

	    QString *local =
	      mp_polymer->monomerText(coordinates->start(),
				       coordinates->end(),
				       true);
	    text->append(*local);

	    // If this is a cross-linked oligomer and we are not appending
	    // text for the __last__ oligomer, then append "~" to let the
	    // user know that the sequences are cross-linked.
	    if (oligomerCount > 1 && iter < oligomerCount - 1)
	      text->append("~");

	    delete(local);
	  }
      }
    else
      {
	if(m_monomerText.size())
	  {
	    *text = m_monomerText;

	    return text;
	  }
	else
	  {
	    if (m_monomerList.size())
	      return Sequence::monomerText(0, m_monomerList.size(), true);
	  }
      }

    return text;
  }


  //! Calculates the masses(mono and avg).
  /*! The calculation is performed by computing the mono and avg masses
    of the sequence stretch as described by the start and end indices in
    the polymer sequence. Default calculation options are used.

    \return true if calculations were successful, false otherwise.

    \sa calculateMasses(const CalcOptions &calcOptions, const
    IonizeRule &ionizeRule)

    \sa Polymer::ionize()

    \sa CalcOptions::CalcOptions()
  */
  bool
  Oligomer::calculateMasses()
  {
    CalcOptions calcOptions;

    calcOptions.setCapping(MXT_CAP_NONE);

    IonizeRule rule;

    return calculateMasses(&calcOptions, &rule);
  }


  //! Calculates the masses(mono and avg).
  /*! The calculation is performed by computing the mono and avg masses
    of the sequence stretch as described by the start and end indices in
    the polymer sequence.

    \param calcOptions Calculation options to be used for the mass
    calculation.

    \param ionizeRule Ionization rule to be used for the mass
    calculation.

    \return true if calculations were successful, false otherwise.

    \sa Polymer::calculateMasses()

    \sa Polymer::ionize()
  */
  bool
  Oligomer::calculateMasses(const CalcOptions *calcOptions,
			     const IonizeRule *ionizeRule)
  {
    Q_ASSERT(calcOptions);

    CalcOptions localOptions(*calcOptions);

    // The coordinates of the oligomer are the following:

    // MAMISGMSGRKAS

    // For a tryptic peptide obtained from protein above, we'd have

    // MAMISGMSGR, that is in the oligomer coordinates:

    // [0] MAMISGMSGR [9]

    // When computing the mass of the oligomer, we have to do a

    // for (iter == [0] ; iter < [9 + 1]; ++iter)

    // Which is why we increment add 1 to m_endIndex in the function below.

    // A polymer might be something made of more than one residual chain
    // in case it is a cross-linked oligomer. Compute the mass fore each
    // residual chain, without accounting for the cross-links...

    m_mono = 0;
    m_avg = 0;

    // An oligomer _is_ a CoordinateList, and we need that list in the
    // calcOptions so that we can call the Polymer::accountMasses
    // function.

    localOptions.setCoordinateList(*this);

    // We do not want to take into account the cross-links because
    // we'll be doing this here and because it cannot work if the
    // cross-links are taken into account from the polymer.

    int flags = localOptions.monomerEntities();
    flags &= ~MXT_MONOMER_CHEMENT_CROSS_LINK;
    localOptions.setMonomerEntities(flags);

    Polymer::accountMasses(mp_polymer,
			    localOptions, &m_mono, &m_avg);

//     qDebug() << __FILE__ << __LINE__
// 	      << "After accounting masses(prior to cross-links):"
// 	      << "mono mass is:" << m_mono;

    // At this point, we have added the mass of each constituent
    // oligomer's residual chain. Let's deal with the cross-links.

    for (int iter = 0; iter < m_crossLinkList.size(); ++iter)
      {
	CrossLink *crossLink = m_crossLinkList.at(iter);

	if(!crossLink->accountMasses(&m_mono, &m_avg))
	  return false;

// 	qDebug() << __FILE__ << __LINE__
// 		  << "After accounting cross-link:"
// 		  << crossLink->name()
// 		  << "mono mass is:" << m_mono;
      }

    // If an ionizeRule is provided, use it. Otherwise, ionize
    // automatically using the ::Ionizable IonizeRule.
    if (ionizeRule)
      {
	// Note that if ionizeRule is invalid, then the ionization is
	// not performed.

	if(ionizeRule->isValid())
	  {
	    /*
	      if (ionize(mp_polymer, *ionizeRule) == -1)
	      The line above is a huge bug. While we should be
	      ionizing this oligomer, we end up ionizing the polymer !
	    */

	    if (ionize(*ionizeRule) == -1)
	      return false;
	  }
      }
    else
      {
	if(ionize() == -1)
	  return false;
      }

//     qDebug() << __FILE__ << __LINE__
// 	      << "Coming out from the calculateMasses function:"
// 	      << "mono mass is:" << m_mono;

    return true;
  }


  //! Tells if the oligomer contains at least one modified monomer.
  /*! Iterates in the polymer sequence from start index to end index and
    if one of the iterated monomers contain a PxmModif modification,
    returns true;

    \return true if at least one monomer is modified, false otherwise.
  */
  bool
  Oligomer::isModified()
  {
    int oligomerCount = CoordinateList::size();

    for (int iter = 0; iter < oligomerCount; ++iter)
      {
	Coordinates *coordinates = CoordinateList::at(iter);

	for(int jter = coordinates->start();
	     jter < coordinates->end() + 1; ++jter)
	  {
// 	    qDebug() << __FILE__ << __LINE__ << "Going to call at() with value"
// 		      << iter;

	    const Monomer *monomer = mp_polymer->at(jter);

	    if (monomer->isModified())
	      return true;
	  }
      }

    return false;
  }


  int
  Oligomer::size()
  {
    int sum = 0;

    int oligomerCount = CoordinateList::size();

    // The size of an oligomer is the sum of all its oligomeric
    // components as described by the various coordinates.

    for (int iter = 0; iter < oligomerCount; ++iter)
      {
	Coordinates *coordinates = CoordinateList::at(iter);

	sum += coordinates->length();
      }

    return sum;
  }


  bool
  Oligomer::encompasses(int index) const
  {
    int oligomerCount = CoordinateList::size();

    for (int iter = 0; iter < oligomerCount; ++iter)
      {
	Coordinates *coordinates = CoordinateList::at(iter);

	if(index <= coordinates->start() && index >= coordinates->end())
	  return true;
      }

    return false;
  }


  bool
  Oligomer::encompasses(const Monomer *monomer) const
  {
    int oligomerCount = CoordinateList::size();

    for (int iter = 0; iter < oligomerCount; ++iter)
      {
	Coordinates *coordinates = CoordinateList::at(iter);

	for(int jter = coordinates->start();
	     jter < coordinates->end() + 1; ++jter)
	  {
// 	    qDebug() << __FILE__ << __LINE__ << "Going to call at() with value"
// 		      << jter;

	    if (mp_polymer->at(jter) == monomer)
	      return true;
	  }
      }

    return false;
  }

} // namespace massXpert