text: /*  $Id: cigar_formatter.cpp 562945 2018-04-30 16:39:32Z kans $
  * ===========================================================================
  *
  *                            PUBLIC DOMAIN NOTICE
  *               National Center for Biotechnology Information
  *
  *  This software/database is a "United States Government Work" under the
  *  terms of the United States Copyright Act.  It was written as part of
  *  the author's official duties as a United States Government employee and
  *  thus cannot be copyrighted.  This software/database is freely available
  *  to the public for use. The National Library of Medicine and the U.S.
  *  Government have not placed any restriction on its use or reproduction.
  *
  *  Although all reasonable efforts have been taken to ensure the accuracy
  *  and reliability of the software and data, the NLM and the U.S.
  *  Government do not and cannot warrant the performance or results that
  *  may be obtained by using this software or data. The NLM and the U.S.
  *  Government disclaim all warranties, express or implied, including
  *  warranties of performance, merchantability or fitness for any particular
  *  purpose.
  *
  *  Please cite the author in any work or product based on this material.
  *
  * ===========================================================================
  *
  * Authors:  Aaron Ucko, Aleksey Grichenko
  *
  * File Description:
  *   Base class for CIGAR formatters.
  *
  * ===========================================================================
  */
  
  #include <ncbi_pch.hpp>
  #include <objects/seqalign/Dense_seg.hpp>
  #include <objects/seqalign/Spliced_seg.hpp>
  #include <objects/seqalign/Seq_align_set.hpp>
  #include <objects/seqalign/Score.hpp>
  #include <objmgr/util/sequence.hpp>
  #include <objtools/error_codes.hpp>
  #include <objtools/format/flat_expt.hpp>
  #include <objtools/format/cigar_formatter.hpp>
  
  
  #define NCBI_USE_ERRCODE_X   Objtools_Fmt_CIGAR
  
  BEGIN_NCBI_SCOPE
  BEGIN_SCOPE(objects)
  
  
  CCIGAR_Formatter::CCIGAR_Formatter(const CSeq_align&    aln,
                                     CScope*              scope,
                                     TCIGARFlags          flags)
      : m_Align(aln),
        m_CurAlign(NULL),
        m_Scope(scope),
        m_Flags(flags),
        m_IsFirstSubalign(true),
        m_IsTrivial(true),
        m_LastType(0),
        m_Frame(-1),
        m_RefRow(-1),
        m_RefSign(1),
        m_TargetRow(-1),
        m_TargetSign(1),
        m_FormatBy(eFormatBy_NotSet)
  {
  }
  
  
  CCIGAR_Formatter::~CCIGAR_Formatter(void)
  {
  }
  
  
  void CCIGAR_Formatter::x_FormatAlignmentRows(void)
  {
      StartAlignment();
      x_FormatAlignmentRows(GetSeq_align(), false);
      EndAlignment();
  }
  
  
  void CCIGAR_Formatter::x_FormatAlignmentRows(const CSeq_align& sa,
                                               bool              width_inverted)
  {
      switch (sa.GetSegs().Which()) {
      case CSeq_align::TSegs::e_Denseg:
          x_FormatDensegRows(sa.GetSegs().GetDenseg(), width_inverted);
          break;
  
      case CSeq_align::TSegs::e_Spliced:
          {
              CRef<CSeq_align> sa2;
              try {
                  sa2 = sa.GetSegs().GetSpliced().AsDiscSeg();
                  if (sa.IsSetScore()) {
                      sa2->SetScore().insert(sa2->SetScore().end(),
                          sa.GetScore().begin(),
                          sa.GetScore().end());
                  }
              } STD_CATCH_ALL_X(1, "CCIGAR_Formatter::x_FormatAlignmentRows")
              if (sa2) {
                  // Conversion from Spliced to Disc inverts meaning of width!!!
                  x_FormatAlignmentRows(*sa2, true);
              }
              break;
          }
  
      case CSeq_align::TSegs::e_Std:
          {
              CRef<CSeq_align> sa2;
              try {
                  sa2 = sa.CreateDensegFromStdseg();
              } STD_CATCH_ALL_X(1, "CCIGAR_Formatter::x_FormatAlignmentRows")
              if (sa2.NotEmpty()  &&  sa2->GetSegs().IsDenseg()) {
                  x_FormatDensegRows(sa2->GetSegs().GetDenseg(), width_inverted);
              }
              break;
          }
  
      case CSeq_align::TSegs::e_Disc:
          {
              ITERATE (CSeq_align_set::Tdata, it, sa.GetSegs().GetDisc().Get()) {
                  m_CurAlign = (*it).GetPointer();
                  StartSubAlignment();
                  x_FormatAlignmentRows(**it, width_inverted);
                  EndSubAlignment();
                  m_CurAlign = NULL;
                  m_IsFirstSubalign = false;
              }
              break;
          }
  
      default: // dendiag or packed; unsupported
          NCBI_THROW(CFlatException, eNotSupported,
              "Conversion of alignments of type dendiag and packed "
              "not supported in current CIGAR output");
      }
  }
  
  
  CCIGAR_Formatter::TNumrow CCIGAR_Formatter::x_GetRowById(const CSeq_id& id)
  {
      CScope* scope = GetScope();
      for (TNumrow row = 0; row < m_AlnMap->GetNumRows(); ++row) {
          if (sequence::IsSameBioseq(m_AlnMap->GetSeqId(row), id, scope)) {
              return row;
          }
      }
      ERR_POST_X(1, "CCIGAR_Formatter::x_GetRowById: "
          "no row with a matching ID found: " << id.AsFastaString());
      return -1;
  }
  
  
  void CCIGAR_Formatter::x_FormatLine(bool width_inverted)
  {
      if (m_TargetRow == m_RefRow) {
          return;
      }
      CNcbiOstrstream cigar;
      m_LastType = 0;
      TSeqPos last_count = 0;
  
      if ( !m_RefId ) {
          m_RefId.Reset(&m_AlnMap->GetSeqId(m_RefRow));
          AdjustSeqIdType(m_RefId);
      }
      if ( !m_TargetId ) {
          m_TargetId.Reset(&m_AlnMap->GetSeqId(m_TargetRow));
          AdjustSeqIdType(m_TargetId);
      }
  
      m_RefWidth =
          (static_cast<size_t>(m_RefRow) < m_DenseSeg->GetWidths().size()) ?
          m_DenseSeg->GetWidths()[m_RefRow] : 1;
      m_RefSign = m_AlnMap->StrandSign(m_RefRow);
      m_TargetWidth =
          (static_cast<size_t>(m_TargetRow) < m_DenseSeg->GetWidths().size()) ?
          m_DenseSeg->GetWidths()[m_TargetRow] : 1;
      m_TargetSign = m_AlnMap->StrandSign(m_TargetRow);
      m_IsTrivial = true;
      TSignedSeqPos last_frameshift = 0;
  
      if (! width_inverted  &&  (m_RefWidth != 1  ||  m_TargetWidth != 1)) {
          // Supporting widths ONLY in the unamiguous case when we
          // know they are WRONG and put there incorrectly from conversion
          // from Spliced-seg. If we didn't get widths that way, we don't
          // know what they mean, so punt if not all widths are 1.
          NCBI_THROW(CFlatException, eNotSupported,
              "Widths in alignments do not have clear semantics, "
              "and thus are not supported in current CIGAR output");
      }
  
      StartRow();
  
      // Is the following correct???
      //
      // Expecting all coordinates to be normalized relative to
      // some reference width, which might be the
      // Least Common Multiple of length in nucleotide bases of
      // the coordinate system used for each row, e.g. using
      // LCM of 3 if either row is protein. The Least Common Multiple
      // would allow accurately representing frameshifts in either
      // sequence.
      //
      // What does width for an alignment really mean, biologically?
      // It can't have arbitrary meaning, because CIGAR has fixed
      // semantics that M/I/D/F/R are in 3-bp units (i.e. one aa) for
      // proteins and 1-bp units for cDNA.
      //
      // Thus, in practice, I think we are expecting widths to be
      // one of (1, 1) for nuc-nuc, (1, 3) for nuc-prot,
      // (3, 1) for prot-nuc, and (3, 3) for prot-prot.
      TSeqPos width = max(m_RefWidth, m_TargetWidth);
      CAlnMap::TNumchunk numseg = m_AlnMap->GetNumSegs();
      for (CAlnMap::TNumchunk idx = 0; idx < numseg; ++idx) {
          // If ref is on minus strand, revert the order of segments.
          CAlnMap::TNumchunk seg = m_RefSign > 0 ? idx : numseg - idx - 1;
          TRange ref_piece = m_AlnMap->GetRange(m_RefRow, seg);
          TRange tgt_piece = m_AlnMap->GetRange(m_TargetRow, seg);
          CAlnMap::TSegTypeFlags ref_flags = m_AlnMap->GetSegType(m_RefRow, seg);
          CAlnMap::TSegTypeFlags tgt_flags = m_AlnMap->GetSegType(m_TargetRow, seg);
          //The type and count are guaranteed set by one of the if/else cases below.  
          char type = 'X'; // Guaranteed set. Pacify compiler.
          TSeqPos count = 0;  // Guaranteed set. Pacify compiler.
          TSignedSeqPos frameshift = 0;
  
          if ( (tgt_flags & CAlnMap::fSeq)  &&
              !(ref_flags & CAlnMap::fSeq) ) {
              // TODO: Handle non-initial protein gap that does not start
              //       on an aa boundary.
              //
              type = 'I';
              if (idx == 0  &&  IsSetFlag(fCIGAR_GffForFlybase)  &&  m_TargetWidth == 3) {
                  // See comments about frame and phase, below.
                  m_Frame = tgt_piece.GetFrom() % m_TargetWidth;
              }
              count = tgt_piece.GetLength()/width;
              frameshift = -(tgt_piece.GetLength()%TSignedSeqPos(width));
              tgt_piece.SetFrom(tgt_piece.GetFrom()/m_TargetWidth);
              tgt_piece.SetTo(tgt_piece.GetTo()/m_TargetWidth);
              m_TargetRange += tgt_piece;
          }
          else if (! (tgt_flags & CAlnMap::fSeq)  &&
              (ref_flags & CAlnMap::fSeq)) {
              // TODO: Handle gap that does not start on an aa boundary.
              //
              type = 'D';
              if (idx == 0  &&  IsSetFlag(fCIGAR_GffForFlybase)  &&  m_RefWidth == 3) {
                  // See comments about frame and phase, below.
                  m_Frame = ref_piece.GetFrom() % m_RefWidth;
              }
              count = ref_piece.GetLength()/width;
              frameshift = +(ref_piece.GetLength()%width);
              // Adjusting for start position, converting to natural cordinates
              // (aa for protein locations, which would imply divide by 3).
              ref_piece.SetFrom(ref_piece.GetFrom()/m_RefWidth);
              ref_piece.SetTo(ref_piece.GetTo()/m_RefWidth);
              m_RefRange += ref_piece;
          }
          else if ((tgt_flags & CAlnMap::fSeq)  &&
              (ref_flags & CAlnMap::fSeq)) {
              // Hanlde case when sequences aligned.
              // The remaining case is when both don't align at all,
              // which shouldn't happen in a pairwise alignment. If we
              // happen to have a multiple alignment, the remaining case
              // would be one that aligns unrelated sequences, thus has
              // no affect on the current GFF3 output.
              // TODO: Resolve why the following implementation is different
              //       from the above historic implementation. The difference
              //       will be in rounding down vs up on single or last
              //       segment.
              //
              type = 'M';
              if (ref_piece.GetLength()  !=  tgt_piece.GetLength()) {
                  // There's a frameshift.. somewhere. Is this valid? Bail.
                  NCBI_THROW(CFlatException, eNotSupported,
                      "Frameshift(s) in Spliced-exon-chunk's diag "
                      "not supported in current CIGAR output");
              }
              if (idx == 0  &&  IsSetFlag(fCIGAR_GffForFlybase)) {
                  // Semantics of the phase aren't defined in GFF3 for
                  // feature types other than a CDS, and this is an alignment.
                  //
                  // Since phase is not required for alignment features, don't
                  // emit one, unless we have been requested with the special
                  // Flybase variant of GFF3 -- they did ask for phase.
                  //
                  // Also, phase can only be interpreted if we have an alignment
                  // in terms of protein aa, and a width of 3 for one or
                  // the other.
                  //
                  // For an alignment, the meaning of phase is ambiguous,
                  // particularly in dealing with a protein-protein
                  // alignment (if ever it allowed alignment to parts of
                  // a codon), and when the seqid is the protein, rather
                  // than the target.
                  //
                  // A protein won't be "reverse complemented" thus,
                  // can assume that it's plus-strand and look at start
                  // position.
                  //
                  // The computation below is actually for the frame.
                  // The phase is not the same, and will be derived from
                  // the frame.
                  if (m_RefWidth == 3) {
                      m_Frame = ref_piece.GetFrom() % m_RefWidth;
                  } else if (m_TargetWidth == 3) {
                      m_Frame = tgt_piece.GetFrom() % m_TargetWidth;
                  }
              }
              // Adjusting for start position, converting to natural cordinates
              // (aa for protein locations, which would imply divide by 3).
              count = ref_piece.GetLength()/width;
              ref_piece.SetFrom(ref_piece.GetFrom()/m_RefWidth);
              ref_piece.SetTo(ref_piece.GetTo()/m_RefWidth);
              m_RefRange += ref_piece;
              tgt_piece.SetFrom(tgt_piece.GetFrom()/m_TargetWidth);
              tgt_piece.SetTo(tgt_piece.GetTo()/m_TargetWidth);
              m_TargetRange += tgt_piece;
          }
          if (type == m_LastType) {
              last_count += count;
              last_frameshift += frameshift;
          } else {
              if (m_LastType) {
                  if (last_count) {
                      m_IsTrivial = false;
                      AddSegment(cigar, m_LastType, last_count);
                  }
                  if (last_frameshift) {
                      m_IsTrivial = false;
                      AddSegment(cigar,
                          (last_frameshift < 0 ? 'F' : 'R'),
                          abs(last_frameshift));
                  }
              }
              m_LastType = type;
              last_count = count;
              last_frameshift = frameshift;
          }
      }
      CNcbiOstrstream aln_out;
  
      AddSegment(cigar, m_LastType, last_count);
      string cigar_string = CNcbiOstrstreamToString(cigar);
  
      AddRow(cigar_string);
  
      EndRow();
  }
  
  
  void CCIGAR_Formatter::x_FormatDensegRows(const CDense_seg& ds,
                                            bool width_inverted)
  {
      m_DenseSeg.Reset(&ds);
  
      // Frame, as a value of 0, 1, 2, or -1 for undefined.
      // This is NOT the same frame as the frame in ASN.1!
      m_Frame = -1;
  
      // WORKAROUND
      // I do believe there is lack of agreement on what
      // the "widths" of a dense-seg mean -- as multiplier, or as divisor.
      //
      // In CSpliced_seg::s_ExonToDenseg, the widths act as divisors,
      // i.e. width 3 means every 3 units (na) in the alignment
      // correspond to 1 unit (aa) on the sequence.
      //
      // In CAlnMix as witnessed by e.g. x_ExtendDSWithWidths,
      // or even CDense_seg::Validate, the widths are a multiplier,
      // i.e. width 3 means every 1 unit (aa) in the alignment
      // corresponds to an alignment of 3 units (na) on the sequence.
      //
      // These definitions are incompatible.
      // The problem with the latter definition as a multiplier,
      // is that the smallest unit of alignment (in a protein-to-nucleotide)
      // is 1 aa = 3 bp... no opportunity for a frameshift. :-(
      //
      // To compensate (or rather, avoid double-compentating), avoid use
      // of widths, and copy to a temporary alignment, storing the old widths
      // for lookup, but reset them in the temporary alignment.
      const CDense_seg* ds_for_alnmix(&ds);
      CDense_seg ds_no_widths;
      if (width_inverted) {
          ds_no_widths.Assign(ds);
          ds_no_widths.ResetWidths();
          ds_for_alnmix = &ds_no_widths;
      }
  
      m_AlnMap.Reset(new CAlnMap(*ds_for_alnmix));
  
      switch ( m_FormatBy ) {
      case eFormatBy_ReferenceId:
          {
              bool by_id = m_RefId.NotNull();
              if ( by_id ) {
                  m_RefRow = x_GetRowById(*m_RefId);
              }
              else {
                  _ASSERT(m_RefRow >= 0);
                  m_RefId.Reset(&m_AlnMap->GetSeqId(m_RefRow));
              }
              StartRows();
              for (m_TargetRow = 0; m_TargetRow < m_AlnMap->GetNumRows(); ++m_TargetRow) {
                  x_FormatLine(width_inverted);
                  m_TargetId.Reset();
              }
              m_TargetRow = -1;
              if ( by_id ) {
                  m_RefRow = -1;
              }
              else {
                  m_RefId.Reset();
              }
              break;
          }
      case eFormatBy_TargetId:
          {
              bool by_id = m_TargetId.NotNull();
              if ( by_id ) {
                  m_TargetRow = x_GetRowById(*m_TargetId);
              }
              else {
                  _ASSERT(m_TargetRow >= 0);
                  m_TargetId.Reset(&m_AlnMap->GetSeqId(m_TargetRow));
                  AdjustSeqIdType(m_TargetId);
              }
              StartRows();
              for (m_RefRow = 0; m_RefRow < m_AlnMap->GetNumRows(); ++m_RefRow) {
                  x_FormatLine(width_inverted);
                  m_RefId.Reset();
              }
              m_RefRow = -1;
              if ( by_id ) {
                  m_TargetRow = -1;
              }
              else {
                  m_TargetId.Reset();
              }
              break;
          }
      default:
          break;
      }
  
      EndRows();
  
      // Reset all values which have no sence anymore
      m_DenseSeg.Reset();
      m_AlnMap.Reset();
      m_IsTrivial = true;
      m_LastType = 0;
      m_Frame = -1;
      m_RefRange = TRange::GetEmpty();
      m_RefSign = 1;
      m_TargetRange = TRange::GetEmpty();
      m_TargetSign = 1;
  }
  
  
  void CCIGAR_Formatter::FormatByReferenceId(const CSeq_id& ref_id)
  {
      m_FormatBy = eFormatBy_ReferenceId;
      m_RefId.Reset(&ref_id);
      m_TargetId.Reset();
      m_RefRow = -1;
      m_TargetRow = -1;
      x_FormatAlignmentRows();
  }
  
  
  void CCIGAR_Formatter::FormatByTargetId(const CSeq_id& target_id)
  {
      m_FormatBy = eFormatBy_TargetId;
      m_RefId.Reset();
      m_TargetId.Reset(&target_id);
      m_RefRow = -1;
      m_TargetRow = -1;
      x_FormatAlignmentRows();
  }
  
  
  void CCIGAR_Formatter::FormatByReferenceRow(TNumrow ref_row)
  {
      m_FormatBy = eFormatBy_ReferenceId;
      m_RefId.Reset();
      m_TargetId.Reset();
      m_RefRow = ref_row;
      m_TargetRow = -1;
      x_FormatAlignmentRows();
  }
  
  
  void CCIGAR_Formatter::FormatByTargetRow(TNumrow target_row)
  {
      m_FormatBy = eFormatBy_TargetId;
      m_RefId.Reset();
      m_TargetId.Reset();
      m_RefRow = -1;
      m_TargetRow = target_row;
      x_FormatAlignmentRows();
  }
  
  
  void CCIGAR_Formatter::AddSegment(CNcbiOstream& cigar,
                                    char seg_type,
                                    TSeqPos seg_len)
  {
      cigar << seg_len << seg_type;
  }
  
  
  END_SCOPE(objects)
  END_NCBI_SCOPE
  