objects/seqalign/Seq_align.cpp

/* $Id: Seq_align.cpp 618884 2020-10-27 16:28:50Z gouriano $
 * ===========================================================================
 *
 *                            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.
 *
 * ===========================================================================
 *
 * Author:  .......
 *
 * File Description:
 *   .......
 *
 * Remark:
 *   This code was originally generated by application DATATOOL
 *   using specifications from the data definition file
 *   'seqalign.asn'.
 */

// standard includes
#include <ncbi_pch.hpp>
#include <objects/seqalign/seqalign_exception.hpp>
#include <objects/seqalign/Dense_seg.hpp>
#include <objects/seqalign/Dense_diag.hpp>
#include <objects/seqalign/Std_seg.hpp>
#include <objects/seqalign/Spliced_seg.hpp>
#include <objects/seqalign/Spliced_exon.hpp>
#include <objects/seqalign/Spliced_exon_chunk.hpp>
#include <objects/seqalign/Sparse_seg.hpp>
#include <objects/seqalign/Seq_align_set.hpp>
#include <objects/seqalign/Score.hpp>
#include <objects/general/Object_id.hpp>
#include <objects/general/User_object.hpp>
#include <objects/seqloc/Seq_loc.hpp>
#include <objects/seqloc/Seq_interval.hpp>
#include <objects/seqloc/Seq_point.hpp>
#include <objects/seq/seq_loc_mapper_base.hpp>
#include <serial/iterator.hpp>
#include <util/static_set.hpp>

// generated includes
#include <objects/seqalign/Seq_align.hpp>

// generated classes

BEGIN_NCBI_SCOPE

BEGIN_objects_SCOPE // namespace ncbi::objects::

// destructor
CSeq_align::~CSeq_align(void)
{
}


CSeq_align::TDim CSeq_align::CheckNumRows(void) const
{
    switch (GetSegs().Which()) {
    case C_Segs::e_Denseg:
        return GetSegs().GetDenseg().CheckNumRows();

        //case C_Segs::e_Packed:
        //return GetSegs().GetPacked().CheckNumRows();

    case C_Segs::e_Sparse:
        return GetSegs().GetSparse().CheckNumRows();

    case C_Segs::e_Spliced:
        {{
            // spliced seg always has two rows: genomic and protein/transcript
            return 2;
        }}

    case C_Segs::e_Disc:
        {{
            TDim numrows = 0;
            ITERATE (C_Segs::TDisc::Tdata, std_i, GetSegs().GetDisc().Get()) {
                TDim seg_numrows = (*std_i)->CheckNumRows();
                if (numrows) {
                    if ( seg_numrows != numrows ) {
                        NCBI_THROW(CSeqalignException, eInvalidAlignment,
                                   "CSeq_align::CheckNumRows(): Number of rows "
                                   "is not the same for each disc seg.");
                    }
                } else {
                    numrows = seg_numrows;
                }
            }
            return numrows;
        }}

    case C_Segs::e_Std:
        {{
            TDim numrows = 0;
            ITERATE (C_Segs::TStd, std_i, GetSegs().GetStd()) {
                TDim seg_numrows = (*std_i)->CheckNumRows();
                if (numrows) {
                    if (seg_numrows != numrows) {
                        NCBI_THROW(CSeqalignException, eInvalidAlignment,
                                   "CSeq_align::CheckNumRows(): Number of rows "
                                   "is not the same for each std seg.");
                    }
                } else {
                    numrows = seg_numrows;
                }
            }
            return numrows;
        }}

    case C_Segs::e_Dendiag:
        {
            TDim numrows = 0;
            ITERATE (C_Segs::TDendiag, it, GetSegs().GetDendiag()) {
                TDim seg_numrows = (*it)->CheckNumRows();
                if (numrows) {
                    if ( seg_numrows != numrows ) {
                        NCBI_THROW(CSeqalignException, eInvalidAlignment,
                                   "CSeq_align::CheckNumRows(): Number of rows "
                                   "is not the same for each dendiag seg.");
                    }
                } else {
                    numrows = seg_numrows;
                }
            }
            return numrows;
        }

    default:
        NCBI_THROW(CSeqalignException, eUnsupported,
                   "CSeq_align::CheckNumRows() currently does not handle "
                   "this type of alignment");
    }
}


CRange<TSeqPos> CSeq_align::GetSeqRange(TDim row) const
{
    switch (GetSegs().Which()) {
    case C_Segs::e_Denseg:
        return GetSegs().GetDenseg().GetSeqRange(row);

    case C_Segs::e_Sparse:
        return GetSegs().GetSparse().GetSeqRange(row);

    case C_Segs::e_Dendiag:
        {
            CRange<TSeqPos> rng;
            ITERATE (C_Segs::TDendiag, dendiag_i, GetSegs().GetDendiag()) {
                rng.CombineWith((*dendiag_i)->GetSeqRange(row));
            }
            return rng;
        }

    case C_Segs::e_Std:
        {
            CRange<TSeqPos> rng;
            CSeq_id seq_id;
            size_t seg_i = 0;
            ITERATE (C_Segs::TStd, std_i, GetSegs().GetStd()) {
                TDim row_i = 0;
                ITERATE (CStd_seg::TLoc, i, (*std_i)->GetLoc()) {
                    const CSeq_loc& loc = **i;
                    if (row_i++ == row) {
                        if (loc.IsInt()) {
                            if ( !seg_i ) {
                                seq_id.Assign(loc.GetInt().GetId());
                            } else if (seq_id.Compare(loc.GetInt().GetId())
                                       != CSeq_id::e_YES) {
                                NCBI_THROW(CSeqalignException,
                                           eInvalidRowNumber,
                                           "CSeq_align::GetSeqRange():"
                                           " Row seqids not consistent."
                                           " Cannot determine range.");
                            }
                            rng.CombineWith
                                (CRange<TSeqPos>
                                 (loc.GetInt().GetFrom(),
                                  loc.GetInt().GetTo()));
                        } else if (loc.IsPnt()) {
                            if ( !seg_i ) {
                                seq_id.Assign(loc.GetPnt().GetId());
                            } else if (seq_id.Compare(loc.GetPnt().GetId())
                                       != CSeq_id::e_YES) {
                                NCBI_THROW(CSeqalignException,
                                           eInvalidRowNumber,
                                           "CSeq_align::GetSeqRange():"
                                           " Row seqids not consistent."
                                           " Cannot determine range.");
                            }
                            rng.CombineWith
                                (CRange<TSeqPos>
                                 (loc.GetPnt().GetPoint(),
                                  loc.GetPnt().GetPoint()));
                        }
                    }
                }
                if (row < 0  ||  row >= row_i) {
                    NCBI_THROW(CSeqalignException, eInvalidRowNumber,
                               "CSeq_align::GetSeqRange():"
                               " Invalid row number");
                }
                if (CanGetDim()  &&  row_i != GetDim()) {
                    NCBI_THROW(CSeqalignException, eInvalidAlignment,
                               "CSeq_align::GetSeqRange():"
                               " loc.size is inconsistent with dim");
                }
                seg_i++;
            }
            if (rng.Empty()) {
                NCBI_THROW(CSeqalignException, eInvalidAlignment,
                           "CSeq_align::GetSeqRange(): Row is empty");
            }
            return rng;
        }

    case C_Segs::e_Disc:
        {
            CRange<TSeqPos> rng;
            ITERATE (C_Segs::TDisc::Tdata, disc_i, GetSegs().GetDisc().Get()) {
                rng.CombineWith((*disc_i)->GetSeqRange(row));
            }
            return rng;
        }

    case C_Segs::e_Spliced:
        return GetSegs().GetSpliced().GetSeqRange(row);

    default:
        NCBI_THROW(CSeqalignException, eUnsupported,
                   "CSeq_align::GetSeqRange() currently does not handle "
                   "this type of alignment.");
    }
}

TSeqPos CSeq_align::GetSeqStart(TDim row) const
{
    switch (GetSegs().Which()) {
    case C_Segs::e_Denseg:
        return GetSegs().GetDenseg().GetSeqStart(row);
    case C_Segs::e_Sparse:
        return GetSegs().GetSparse().GetSeqStart(row);
    case C_Segs::e_Spliced:
        return GetSegs().GetSpliced().GetSeqStart(row);
    case C_Segs::e_Dendiag:
    case C_Segs::e_Std:
    case C_Segs::e_Disc:
        return GetSeqRange(row).GetFrom();
    default:
        NCBI_THROW(CSeqalignException, eUnsupported,
                   "CSeq_align::GetSeqStart() currently does not handle "
                   "this type of alignment.");
    }
}


TSeqPos CSeq_align::GetSeqStop (TDim row) const
{
    switch (GetSegs().Which()) {
    case C_Segs::e_Denseg:
        return GetSegs().GetDenseg().GetSeqStop(row);
    case C_Segs::e_Sparse:
        return GetSegs().GetSparse().GetSeqStop(row);
    case C_Segs::e_Spliced:
        return GetSegs().GetSpliced().GetSeqStop(row);
    case C_Segs::e_Dendiag:
    case C_Segs::e_Std:
    case C_Segs::e_Disc:
        return GetSeqRange(row).GetTo();
    default:
        NCBI_THROW(CSeqalignException, eUnsupported,
                   "CSeq_align::GetSeqStop() currently does not handle "
                   "this type of alignment.");
    }
}


ENa_strand CSeq_align::GetSeqStrand(TDim row) const
{
    switch (GetSegs().Which()) {
    case C_Segs::e_Denseg:
        return GetSegs().GetDenseg().GetSeqStrand(row);
    case C_Segs::e_Sparse:
        return GetSegs().GetSparse().GetSeqStrand(row);
    case C_Segs::e_Spliced:
        return GetSegs().GetSpliced().GetSeqStrand(row);
    case C_Segs::e_Disc:
        return GetSegs().GetDisc().Get().front()->GetSeqStrand(row);
    case C_Segs::e_Std:
        return GetSegs().GetStd().front()->GetLoc()[row]->GetStrand();
    case C_Segs::e_Dendiag:
        return GetSegs().GetDendiag().front()->GetSeqStrand(row);
    default:
        NCBI_THROW(CSeqalignException, eUnsupported,
                   "CSeq_align::GetSeqStrand() currently does not handle "
                   "this type of alignment.");
    }
}


const CSeq_id& CSeq_align::GetSeq_id(TDim row) const
{
    switch (GetSegs().Which()) {
    case C_Segs::e_Denseg:
        return GetSegs().GetDenseg().GetSeq_id(row);

    case C_Segs::e_Sparse:
        return GetSegs().GetSparse().GetSeq_id(row);

    case C_Segs::e_Dendiag:
        {{
            // If segments have different number of rows, this will try
            // to find the segment with enough rows to get the id.
            ITERATE(CSeq_align::C_Segs::TDendiag, seg, GetSegs().GetDendiag()) {
                if( (*seg)->IsSetIds()  &&
                    (size_t)row < (*seg)->GetIds().size() ) {
                    return *(*seg)->GetIds()[row];
                }
            }
            break;
        }}

    case C_Segs::e_Std:
        {{
            // If segments have different number of rows, this will try
            // to find the segment with enough rows to get the id.
            ITERATE(CSeq_align::C_Segs::TStd, seg, GetSegs().GetStd()) {
                if ( (*seg)->IsSetIds()  &&
                     (size_t)row < (*seg)->GetIds().size()) {
                    return *((*seg)->GetIds()[row]);
                }
                else if ( (*seg)->IsSetLoc()  &&
                     (size_t)row < (*seg)->GetLoc().size() ) {
                    const CSeq_loc& loc = *(*seg)->GetLoc()[row];
                    CConstRef<CSeq_id> id(loc.GetId());
                    if (id) {
                        return *id;
                    }
                }
            }
        }}
        break;

    case C_Segs::e_Disc:
        {{
            // Try to find a sub-alignment for which we can get a
            // Seq-id for this row.
            ITERATE (CSeq_align_set::Tdata, sub_aln, GetSegs().GetDisc().Get()) {
                try {
                    const CSeq_id& rv = (*sub_aln)->GetSeq_id(row);
                    return rv;
                }
                catch (const CSeqalignException&) {
                }
            }
        }}
        break;

    case C_Segs::e_Spliced:
        {{
            // Technically, there is no row order for product and product.
            // However, in spliced seg, since product comes first, we assign it
            // as row 0.
            // Hence, the genomic is assigned to row 1.
            const C_Segs::TSpliced& spliced_seg = GetSegs().GetSpliced();
            if (row == 0 && spliced_seg.IsSetProduct_id()) {
                return spliced_seg.GetProduct_id();
            } else if (row == 1 && spliced_seg.IsSetGenomic_id()) {
                return spliced_seg.GetGenomic_id();
            }
        }}
        break;

    default:
        NCBI_THROW(CSeqalignException, eUnsupported,
                   "CSeq_align::GetSeq_id() currently does not handle "
                   "this type of alignment.");
    }

    NCBI_THROW(CSeqalignException, eInvalidRowNumber,
               "CSeq_align::GetSeq_id(): "
               "can not get seq-id for the row requested.");
    // return CSeq_id();
}


typedef SStaticPair<CSeq_align::EScoreType, const char*> TScoreNamePair;
static const TScoreNamePair sc_ScoreNames[] = {
    { CSeq_align::eScore_Score,           "score" },
    { CSeq_align::eScore_Blast,           "score" },
    { CSeq_align::eScore_BitScore,        "bit_score" },
    { CSeq_align::eScore_EValue,          "e_value" },
    { CSeq_align::eScore_AlignLength,     "align_length" },
    { CSeq_align::eScore_IdentityCount,   "num_ident" },
    { CSeq_align::eScore_PositiveCount,   "num_positives" },
    { CSeq_align::eScore_NegativeCount,   "num_negatives" },
    { CSeq_align::eScore_MismatchCount,   "num_mismatch" },
    { CSeq_align::eScore_GapCount,        "gap_count" },
    { CSeq_align::eScore_PercentIdentity_Gapped, "pct_identity_gap" },
    { CSeq_align::eScore_PercentIdentity_Ungapped, "pct_identity_ungap" },
    { CSeq_align::eScore_PercentIdentity_GapOpeningOnly, "pct_identity_gapopen_only" },
    { CSeq_align::eScore_PercentCoverage, "pct_coverage" },
    { CSeq_align::eScore_SumEValue,       "sum_e" },
    { CSeq_align::eScore_CompAdjMethod,   "comp_adjustment_method" },
    { CSeq_align::eScore_HighQualityPercentCoverage, "pct_coverage_hiqual" },
    { CSeq_align::eScore_Matches,         "matches"},
    { CSeq_align::eScore_OverallIdentity, "identity"},
    { CSeq_align::eScore_Splices,         "splices"},
    { CSeq_align::eScore_ConsensusSplices, "consensus_splices"},
    { CSeq_align::eScore_ProductCoverage, "product_coverage"},
    { CSeq_align::eScore_ExonIdentity,    "exon_identity"}
};

static const char* const sc_ScoreHelpText[] = {
    "Blast score",
    "Blast score",
    "Blast-style bit score",
    "Blast-style e-value",
    "Length of the aligned segments, including the length of all gap segments",
    "Count of identities",
    "Count of positives; protein-to-DNA score",
    "Count of negatives; protein-to-DNA score",
    "Count of mismatches",
    "Number of gaps in the alignment",
    "Percent identity (0.0-100.0); count each base in a gap in any row as a mismatch",
    "Percent identity (0.0-100.0); don't count gaps",
    "Percent identity (0.0-100.0); count a gap of any length in any row as a mismatch of length 1",
    "Percentage of query sequence aligned to subject (0.0-100.0)",
    "Blast-style sum_e",
    "Composition-adjustment method from BLAST",
    "Percent coverage (0.0-100.0) of high quality region",
    "Count of identities",
    "Percent identity; count gaps within exons and introns on product",
    "Number of splices; 2 x number of introns that have no gap on product",
    "Number of splices with consensus splice sequence",
    "Percentage of query sequence aligned to subject (0.0-100.0)",
    "Percent identity; count gaps within exons only"
};

static const bool sc_IsInteger[] = {
    true,   // eScore_Score
    true,   // eScore_Blast
    false,  // eScore_BitScore
    false,  // eScore_EValue
    true,   // eScore_AlignLength
    true,   // eScore_IdentityCount
    true,   // eScore_PositiveCount
    true,   // eScore_NegativeCount
    true,   // eScore_MismatchCount
    true,   // eScore_GapCount
    false,  // eScore_PercentIdentity_Gapped
    false,  // eScore_PercentIdentity_Ungapped
    false,  // eScore_PercentIdentity_GapOpeningOnly
    false,  // eScore_PercentCoverage
    false,  // eScore_SumEValue
    false,  // eScore_CompAdjMethod
    false,  // eScore_HighQualityPercentCoverage
    true,   // eScore_Matches
    false,  // eScore_OverallIdentity
    true,   // eScore_Splices
    true,   // eScore_ConsensusSplices
    false,  // eScore_ProductCoverage
    false   // eScore_ExonIdentity
};


const CSeq_align::TScoreNameMap &CSeq_align::ScoreNameMap()
{
    DEFINE_STATIC_FAST_MUTEX(s_ScoreNameMap_mutex);
    CFastMutexGuard guard(s_ScoreNameMap_mutex);
    {
        static TScoreNameMap m_ScoreNameMap;
        if (m_ScoreNameMap.empty()) {
            /// initialize map
            for(unsigned score = eScore_Blast;
                score <= eScore_ExonIdentity; ++score)
                {
                    m_ScoreNameMap[sc_ScoreNames[score].second] =
                        sc_ScoreNames[score].first;
                }
        }

        return m_ScoreNameMap;
    }
}

string CSeq_align::ScoreName(EScoreType score)
{
    return sc_ScoreNames[score].second;
}

string CSeq_align::HelpText(EScoreType score)
{
    return sc_ScoreHelpText[score];
}

bool CSeq_align::IsIntegerScore(EScoreType score)
{
    return sc_IsInteger[score];
}


/// retrieve a named score object
CConstRef<CScore> CSeq_align::x_GetNamedScore(const string& name) const
{
    CConstRef<CScore> score;
    if (IsSetScore()) {
        ITERATE (TScore, iter, GetScore()) {
            if ((*iter)->IsSetId()  &&  (*iter)->GetId().IsStr()  &&
                (*iter)->GetId().GetStr() == name) {
                score = *iter;
                break;
            }
        }
    }

    return score;
}


CRef<CScore> CSeq_align::x_SetNamedScore(const string& name)
{
    CRef<CScore> score;
    if (IsSetScore()) {
        NON_CONST_ITERATE (TScore, iter, SetScore()) {
            if ((*iter)->IsSetId()  &&  (*iter)->GetId().IsStr()  &&
                (*iter)->GetId().GetStr() == name) {
                score = *iter;
                break;
            }
        }
    }

    if ( !score ) {
        score.Reset(new CScore);
        score->SetId().SetStr(name);
        SetScore().push_back(score);
    }

    return score;
}


///---------------------------------------------------------------------------
/// PRE : name of score to return
/// POST: whether or not we found that score; score converted to int
bool CSeq_align::GetNamedScore(const string &id, int &score) const
{
    CConstRef<CScore> ref = x_GetNamedScore(id);
    if (ref) {
        if (ref->GetValue().IsInt()) {
            score = ref->GetValue().GetInt();
        } else {
            score = (int)ref->GetValue().GetReal();
        }
        return true;
    }
    return false;
}

///---------------------------------------------------------------------------
/// PRE : name of score to return
/// POST: whether or not we found that score; score converted to double
bool CSeq_align::GetNamedScore(const string &id, double &score) const
{
    CConstRef<CScore> ref = x_GetNamedScore(id);
    if (ref) {
        if (ref->GetValue().IsInt()) {
            score = (double)ref->GetValue().GetInt();
        } else {
            score = ref->GetValue().GetReal();
        }
        return true;
    }
    return false;
}


bool CSeq_align::GetNamedScore(EScoreType type, int& score) const
{
    return GetNamedScore(sc_ScoreNames[type].second, score);
}

bool CSeq_align::GetNamedScore(EScoreType type, double& score) const
{
    return GetNamedScore(sc_ScoreNames[type].second, score);
}


void CSeq_align::ResetNamedScore(const string& name)
{
    if (IsSetScore()) {
        NON_CONST_ITERATE (TScore, iter, SetScore()) {
            if ((*iter)->IsSetId()  &&  (*iter)->GetId().IsStr()  &&
                (*iter)->GetId().GetStr() == name) {
                SetScore().erase(iter);
                break;
            }
        }
    }
}

void CSeq_align::ResetNamedScore(EScoreType type)
{
    ResetNamedScore(sc_ScoreNames[type].second);
}

void CSeq_align::SetNamedScore(EScoreType type, int score)
{
    CRef<CScore> ref = x_SetNamedScore(sc_ScoreNames[type].second);
    ref->SetValue().SetInt(score);
}

void CSeq_align::SetNamedScore(EScoreType type, double score)
{
    CRef<CScore> ref = x_SetNamedScore(sc_ScoreNames[type].second);
    ref->SetValue().SetReal(score);
}

void CSeq_align::SetNamedScore(const string& id, int score)
{
    CRef<CScore> ref = x_SetNamedScore(id);
    ref->SetValue().SetInt(score);
}

void CSeq_align::SetNamedScore(const string& id, double score)
{
    CRef<CScore> ref = x_SetNamedScore(id);
    ref->SetValue().SetReal(score);
}


void CSeq_align::Validate(bool full_test) const
{
    switch (GetSegs().Which()) {
    case TSegs::e_Dendiag:
        ITERATE(TSegs::TDendiag, dendiag_it, GetSegs().GetDendiag()) {
            (*dendiag_it)->Validate();
        }
        break;
    case C_Segs::e_Denseg:
        GetSegs().GetDenseg().Validate(full_test);
        break;
    case C_Segs::e_Disc:
        ITERATE(CSeq_align_set::Tdata, seq_align_it, GetSegs().GetDisc().Get()) {
            (*seq_align_it)->Validate(full_test);
        }
        break;
    case C_Segs::e_Std:
        CheckNumRows();
        break;
    case C_Segs::e_Spliced:
        GetSegs().GetSpliced().Validate(full_test);
        break;
    case C_Segs::e_Sparse:
        GetSegs().GetSparse().Validate(full_test);
        break;
    default:
        NCBI_THROW(CSeqalignException, eUnsupported,
                   "CSeq_align::Validate() currently does not handle "
                   "this type of alignment");
    }
}


///---------------------------------------------------------------------------
/// PRE : currently only implemented for dense-seg segments
/// POST: same alignment, opposite orientation
void CSeq_align::Reverse(void)
{
    switch (GetSegs().Which()) {
    case C_Segs::e_Denseg:
        SetSegs().SetDenseg().Reverse();
        break;
    default:
        NCBI_THROW(CSeqalignException, eUnsupported,
                   "CSeq_align::Reverse() currently only handles dense-seg "
                   "alignments");
    }
}

///---------------------------------------------------------------------------
/// PRE : currently only implemented for dense-seg segments; two row numbers
/// POST: same alignment, position of the two rows has been swapped
void CSeq_align::SwapRows(TDim row1, TDim row2)
{
    switch (GetSegs().Which()) {
    case C_Segs::e_Denseg:
        SetSegs().SetDenseg().SwapRows(row1, row2);
        break;

    case C_Segs::e_Std:
        NON_CONST_ITERATE (C_Segs::TStd, it, SetSegs().SetStd()) {
            (*it)->SwapRows(row1, row2);
        }
        break;

    case C_Segs::e_Disc:
        SetSegs().SetDisc().SwapRows(row1, row2);
        break;
    default:
        NCBI_THROW(CSeqalignException, eUnsupported,
                   "CSeq_align::SwapRows currently only handles dense-seg "
                   "alignments");
    }
}

///----------------------------------------------------------------------------
/// PRE : the Seq-align has StdSeg segs
/// POST: Seq_align of type Dense-seg is created with m_Widths if necessary
CRef<CSeq_align>
CSeq_align::CreateDensegFromStdseg(SSeqIdChooser* SeqIdChooser) const
{
    if ( !GetSegs().IsStd() ) {
        NCBI_THROW(CSeqalignException, eInvalidInputAlignment,
                   "CSeq_align::CreateDensegFromStdseg(): "
                   "Input Seq-align should have segs of type StdSeg!");
    }

    CRef<CSeq_align> sa(new CSeq_align);
    sa->SetType(eType_not_set);
    if (IsSetScore()) {
        sa->SetScore() = GetScore();
    }
    CDense_seg& ds = sa->SetSegs().SetDenseg();

    typedef CDense_seg::TDim    TNumrow;
    typedef CDense_seg::TNumseg TNumseg;

    vector<TSeqPos>       row_lens;
    CDense_seg::TLens&    seg_lens = ds.SetLens();
    CDense_seg::TStarts&  starts   = ds.SetStarts();
    CDense_seg::TStrands& strands  = ds.SetStrands();
    CDense_seg::TWidths&  widths   = ds.SetWidths();
    vector<bool>          widths_determined;

    TSeqPos row_len;
    TSeqPos from, to;
    ENa_strand strand;


    TNumseg seg = 0;
    TNumrow dim = 0, row = 0;
    ITERATE (C_Segs::TStd, std_i, GetSegs().GetStd()) {

        const CStd_seg& ss = **std_i;

        seg_lens.push_back(0);
        TSeqPos& seg_len = seg_lens.back();
        row_lens.clear();
        widths_determined.push_back(false);

        row = 0;
        ITERATE (CStd_seg::TLoc, i, ss.GetLoc()) {

            const CSeq_id* seq_id = 0;
            // push back initialization values
            if (seg == 0) {
                widths.push_back(0);
                strands.push_back(eNa_strand_unknown);
            }

            if ((*i)->IsInt()) {
                const CSeq_interval& interval = (*i)->GetInt();

                // determine start and len
                from = interval.GetFrom();
                to = interval.GetTo();
                starts.push_back(from);
                row_len = to - from + 1;
                row_lens.push_back(row_len);

                int width = 0;
                // try to determine/check the seg_len and width
                if (!seg_len) {
                    width = 0;
                    seg_len = row_len;
                } else {
                    if (row_len * 3 == seg_len) {
                        seg_len /= 3;
                        width = 1;
                    } else if (row_len / 3 == seg_len) {
                        width = 3;
                    } else if (row_len != seg_len) {
                        NCBI_THROW(CSeqalignException, eInvalidInputAlignment,
                                   "CreateDensegFromStdseg(): "
                                   "Std-seg segment lengths not accurate!");
                    }
                }
                if (width > 0) {
                    widths_determined[seg] = true;
                    if (widths[row] > 0  &&  widths[row] != width) {
                        NCBI_THROW(CSeqalignException, eInvalidInputAlignment,
                                   "CreateDensegFromStdseg(): "
                                   "Std-seg segment lengths not accurate!");
                    } else {
                        widths[row] = width;
                    }
                }

                // get the id
                seq_id = &(*i)->GetInt().GetId();

                // determine/check the strand
                if (interval.CanGetStrand()) {
                    strand = interval.GetStrand();
                    if (seg == 0  ||  strands[row] == eNa_strand_unknown) {
                        strands[row] = strand;
                    } else {
                        if (strands[row] != strand) {
                            NCBI_THROW(CSeqalignException,
                                       eInvalidInputAlignment,
                                       "CreateDensegFromStdseg(): "
                                       "Inconsistent strands!");
                        }
                    }
                }
            } else if ((*i)->IsEmpty()) {
                starts.push_back(-1);
                if (seg == 0) {
                    strands[row] = eNa_strand_unknown;
                }
                seq_id = &(*i)->GetEmpty();
                row_lens.push_back(0);
            }

            // determine/check the id
            if (seg == 0) {
                CRef<CSeq_id> id(new CSeq_id);
                SerialAssign(*id, *seq_id);
                ds.SetIds().push_back(id);
            } else {
                CSeq_id& id = *ds.SetIds()[row];
                if (!SerialEquals(id, *seq_id)) {
                    if (SeqIdChooser) {
                        SeqIdChooser->ChooseSeqId(id, *seq_id);
                    } else {
                        string errstr =
                            string("CreateDensegFromStdseg(): Seq-ids: ") +
                            id.AsFastaString() + " and " +
                            seq_id->AsFastaString() + " are not identical!" +
                            " Without the OM it cannot be determined if they belong to"
                            " the same sequence."
                            " Define and pass ChooseSeqId to resolve seq-ids.";
                        NCBI_THROW(CSeqalignException, eInvalidInputAlignment, errstr);
                    }
                }
            }

            // next row
            row++;
            if (seg == 0) {
                dim++;
            }
        }
        if (dim != ss.GetDim()  ||  row != dim) {
            NCBI_THROW(CSeqalignException, eInvalidInputAlignment,
                       "CreateDensegFromStdseg(): "
                       "Inconsistent dimentions!");
        }

        if (widths_determined[seg]) {
            // go back and determine/check widths
            for (row = 0; row < dim; row++) {
                if ((row_len = row_lens[row]) > 0) {
                    int width = 0;
                    if (row_len == seg_len * 3) {
                        width = 3;
                    } else if (row_len == seg_len) {
                        width = 1;
                    }
                    if (widths[row] > 0  &&  widths[row] != width) {
                        NCBI_THROW(CSeqalignException, eInvalidInputAlignment,
                                   "CreateDensegFromStdseg(): "
                                   "Std-seg segment lengths not accurate!");
                    } else {
                        widths[row] = width;
                    }
                }
            }
        }

        // next seg
        seg++;
    }

    ds.SetDim(dim);
    ds.SetNumseg(seg);

    // go back and finish lens determination
    bool widths_failure = false;
    bool widths_success = false;
    for (seg = 0; seg < ds.GetNumseg(); seg++) {
        if (!widths_determined[seg]) {
            for(row = 0; row < dim; row++) {
                if (starts[seg * dim + row] >= 0) {
                    int width = widths[row];
                    if (width == 3) {
                        seg_lens[seg] /= 3;
                    } else if (width == 0) {
                        widths_failure = true;
                    }
                    break;
                }
            }
        } else {
            widths_success = true;
        }
    }

    if (widths_failure) {
        if (widths_success) {
            NCBI_THROW(CSeqalignException, eInvalidInputAlignment,
                       "CreateDensegFromStdseg(): "
                       "Some widths cannot be determined!");
        } else {
            ds.ResetWidths();
        }
    }

    // finish the strands
    for (seg = 1; seg < ds.GetNumseg(); seg++) {
        for (row = 0; row < dim; row++) {
            strands.push_back(strands[row]);
        }
    }

    return sa;
}


///----------------------------------------------------------------------------
/// PRE : the Seq-align is a Dense-seg of aligned nucleotide sequences
/// POST: Seq_align of type Dense-seg is created with each of the m_Widths
///       equal to 3 and m_Lengths devided by 3.
CRef<CSeq_align>
CSeq_align::CreateTranslatedDensegFromNADenseg() const
{
    if ( !GetSegs().IsDenseg() ) {
        NCBI_THROW(CSeqalignException, eInvalidInputAlignment,
                   "CSeq_align::CreateTranslatedDensegFromNADenseg(): "
                   "Input Seq-align should have segs of type Dense-seg!");
    }

    CRef<CSeq_align> sa(new CSeq_align);
    sa->SetType(eType_not_set);

    if (GetSegs().GetDenseg().IsSetWidths()) {
        NCBI_THROW(CSeqalignException, eInvalidInputAlignment,
                   "CSeq_align::CreateTranslatedDensegFromNADenseg(): "
                   "Widths already exist for the original alignment");
    }

    // copy from the original
    sa->Assign(*this);

    CDense_seg& ds = sa->SetSegs().SetDenseg();

    // fix the lengths
    const CDense_seg::TLens& orig_lens = GetSegs().GetDenseg().GetLens();
    CDense_seg::TLens&       lens      = ds.SetLens();

    for (CDense_seg::TNumseg numseg = 0; numseg < ds.GetNumseg(); numseg++) {
        if (orig_lens[numseg] % 3) {
            string errstr =
                string("CSeq_align::CreateTranslatedDensegFromNADenseg(): ") +
                "Length of segment " + NStr::IntToString(numseg) +
                " is not divisible by 3.";
            NCBI_THROW(CSeqalignException, eInvalidInputAlignment, errstr);
        } else {
            lens[numseg] = orig_lens[numseg] / 3;
        }
    }

    // add the widths
    ds.SetWidths().resize(ds.GetDim(), 3);

#if _DEBUG
    ds.Validate(true);
#endif

    return sa;
}


/// Strict weak ordering for pairs (by first)
/// Used by CreateDensegFromDisc
template <typename T, typename Pred = less<TSeqPos> >
struct ds_cmp {
    bool operator()(const T& x, const T& y) {
        return m_Pred(x.first, y.first);
    }
private:
    Pred m_Pred;
};


CRef<CSeq_align>
CSeq_align::CreateDensegFromDisc(SSeqIdChooser* SeqIdChooser) const
{
    if ( !GetSegs().IsDisc() ) {
        NCBI_THROW(CSeqalignException, eInvalidInputAlignment,
                   "CSeq_align::CreateDensegFromDisc(): "
                   "Input Seq-align should have segs of type StdSeg!");
    }

    CRef<CSeq_align> new_sa(new CSeq_align);
    new_sa->SetType(eType_not_set);
    if (IsSetScore()) {
        new_sa->SetScore() = GetScore();
    }

    CDense_seg& new_ds = new_sa->SetSegs().SetDenseg();

    new_ds.SetDim(0);
    new_ds.SetNumseg(0);


    /// Order the discontinuous densegs
    typedef pair<TSeqPos, const CDense_seg *> TPosDsPair;
    typedef vector<TPosDsPair> TDsVec;
    TDsVec ds_vec;
    ds_vec.reserve(GetSegs().GetDisc().Get().size());
    int strand = -1;
    ITERATE (CSeq_align_set::Tdata, sa_i, GetSegs().GetDisc().Get()) {
        const CDense_seg& ds = (*sa_i)->GetSegs().GetDenseg();
        ds_vec.push_back(make_pair<TSeqPos, const CDense_seg *>(ds.GetSeqStart(0), &ds));
        if (ds.IsSetStrands()  &&
            !ds.GetStrands().empty()) {
            if (strand < 0) {
                strand = ds.GetStrands()[0];
            } else {
                if (strand != ds.GetStrands()[0]) {
                    NCBI_THROW(CSeqalignException, eInvalidInputAlignment,
                               "CreateDensegFromDisc(): "
                               "Inconsistent strands!");
                }
            }
        }
    }
    if ( !IsReverse(ENa_strand(strand)) ) {
        sort(ds_vec.begin(), ds_vec.end(),
             ds_cmp<TPosDsPair>());
    } else {
        sort(ds_vec.begin(), ds_vec.end(),
             ds_cmp<TPosDsPair, greater<TSeqPos> >());
    }


    /// First pass: determine dim & numseg
    CDense_seg::TStrands single_segment_strands;
    ITERATE(TDsVec, ds_i, ds_vec) {
        const CDense_seg& ds = *ds_i->second;

        /// Numseg
        new_ds.SetNumseg() += ds.GetNumseg();

        /// Dim
        if ( !new_ds.GetDim() ) {
            new_ds.SetDim(ds.GetDim());
        } else if (new_ds.GetDim() != ds.GetDim() ) {
            NCBI_THROW(CSeqalignException, eInvalidInputAlignment,
                       "CreateDensegFromDisc(): "
                       "All disc dense-segs need to have the same dimension!");
        }

        /// Strands?
        if ( !ds.GetStrands().empty() ) {
            if (single_segment_strands.empty()) {
                single_segment_strands.assign(ds.GetStrands().begin(),
                    ds.GetStrands().begin() + ds.GetDim());
            } else {
                if ( !equal(single_segment_strands.begin(),
                            single_segment_strands.end(),
                            ds.GetStrands().begin()) ) {
                    NCBI_THROW(CSeqalignException, eInvalidInputAlignment,
                               "CreateDensegFromDisc(): "
                               "All disc dense-segs need to have the same strands!");
                }
            }
        }
    }

    new_ds.SetStarts().resize(new_ds.GetDim() * new_ds.GetNumseg());
    new_ds.SetLens().resize(new_ds.GetNumseg());
    if ( !single_segment_strands.empty() ) {
        /// Multiply the strands by the number of segments.
        new_ds.SetStrands().reserve(new_ds.GetDim() * new_ds.GetNumseg());
        for (CDense_seg::TNumseg seg = 0;
             seg < new_ds.GetNumseg();
             ++seg) {
            new_ds.SetStrands().insert(new_ds.SetStrands().end(),
                single_segment_strands.begin(),
                single_segment_strands.begin() + new_ds.GetDim());
        }
    }


    /// Second pass: validate and set ids and starts

    CDense_seg::TNumseg new_seg = 0;
    int                 new_starts_i = 0;

    ITERATE(TDsVec, ds_i, ds_vec) {
        const CDense_seg& ds = *ds_i->second;

        _ASSERT(ds.GetStarts().size() == (size_t)ds.GetNumseg() * ds.GetDim());
        _ASSERT(new_ds.GetDim() == ds.GetDim());

        /// Ids
        if (new_ds.GetIds().empty()) {
            new_ds.SetIds().resize(new_ds.GetDim());
            for (CDense_seg::TDim row = 0;  row < ds.GetDim();  ++row) {
                CRef<CSeq_id> id(new CSeq_id);
                SerialAssign(*id, *ds.GetIds()[row]);
                new_ds.SetIds()[row] = id;
            }
        } else {
            if (SeqIdChooser) {
                for (CDense_seg::TDim row = 0;  row < ds.GetDim();  ++row) {
                    SeqIdChooser->ChooseSeqId(*new_ds.SetIds()[row], *ds.GetIds()[row]);
                }
            } else {
#if _DEBUG
                for (CDense_seg::TDim row = 0;  row < ds.GetDim();  ++row) {
                    const CSeq_id& new_id = *new_ds.GetIds()[row];
                    const CSeq_id& id     = *ds.GetIds()[row];
                    if ( !SerialEquals(new_id, id) ) {
                        string errstr =
                            string("CreateDensegFromDisc(): Seq-ids: ") +
                            new_id.AsFastaString() + " and " +
                            id.AsFastaString() + " are not identical!" +
                            " Without the OM it cannot be determined if they belong to"
                            " the same sequence."
                            " Define and pass ChooseSeqId to resolve seq-ids.";
                        NCBI_THROW(CSeqalignException, eInvalidInputAlignment, errstr);
                    }
                }
#endif
            }
        }

        /// Starts
        int starts_i = 0;
        for (CDense_seg::TNumseg seg = 0;
             seg < ds.GetNumseg();
             ++new_seg, ++seg) {

            new_ds.SetLens()[new_seg] = ds.GetLens()[seg];

            for (CDense_seg::TDim row = 0;
                 row < ds.GetDim();
                 ++starts_i, ++new_starts_i, ++row) {
                new_ds.SetStarts()[new_starts_i] = ds.GetStarts()[starts_i];
            }
        }
    }


    /// Perform full test (including segment starts consistency check)
    new_sa->Validate(true);


    return new_sa;
}

static TSeqPos s_Distance(const TSeqRange &range1,
                          const TSeqRange &range2)
{
    if (range1.IntersectingWith(range2)) {
        return 0;
    }
    return max(range1,range2).GetFrom() - min(range1,range2).GetTo();
}

/// Join the alignments in the set into one alignment.
static CRef<CSeq_align> s_GetJoinedAlignment(const CSeq_align_set &aligns)
{
    if (aligns.Get().size() == 1) {
        return aligns.Get().front();
    }

    /// There's more than one, so create a disc alignment
    CRef<CSeq_align> align(new CSeq_align);
    align->SetDim(2);
    align->SetType(CSeq_align::eType_partial);
    align->SetSegs().SetDisc().Set() = aligns.Get();

    try {
        align->Validate(true);
    }
    catch (CException& e) {
        ERR_POST(Error << e);
        cerr << MSerial_AsnText << aligns;
        cerr << MSerial_AsnText << *align;
        throw;
    }
    return align;
}

void CSeq_align::SplitOnLongDiscontinuity(list< CRef<CSeq_align> >& aligns,
                                          TSeqPos discontinuity_threshold) const
{
    if (IsSetDim() && GetDim() > 2) {
        NCBI_THROW(CException, eUnknown,
            "SplitOnLongDiscontinuity() only implemented for pairwise alignments");
    }

    size_t num_splits = 0;
    switch (GetSegs().Which()) {
    case TSegs::e_Denseg:
    {{
        // determine if we can extract any slice based on a discontinuity
        const CDense_seg& seg = GetSegs().GetDenseg();
        CDense_seg::TNumseg last_seg = 0;
        CDense_seg::TNumseg curr_seg = 0;
        TSeqRange curr_range0;
        TSeqRange curr_range1;
        for ( ;  curr_seg < seg.GetNumseg();  ++curr_seg) {
            TSignedSeqPos p0_curr = seg.GetStarts()[ curr_seg * 2 ];
            TSignedSeqPos p1_curr = seg.GetStarts()[ curr_seg * 2 + 1];
            if (p0_curr < 0 || p1_curr < 0) {
                continue;
            }

            TSeqPos curr_seg_len = seg.GetLens()[curr_seg];

            TSeqRange prev_range0 = curr_range0;
            TSeqRange prev_range1 = curr_range1;
            curr_range0 = TSeqRange(p0_curr, p0_curr+curr_seg_len);
            curr_range1 = TSeqRange(p1_curr, p1_curr+curr_seg_len);

            if (prev_range0.Empty()) {
                /// We're at start of alignment, so there's no gap to check for
                continue;
            }

            if (prev_range0.NotEmpty() &&
                (s_Distance(curr_range0, prev_range0) > discontinuity_threshold ||
                 s_Distance(curr_range1, prev_range1) > discontinuity_threshold))
            {
                aligns.push_back(x_CreateSubsegAlignment(last_seg, curr_seg - 1));
                ++num_splits;
                last_seg = curr_seg;
            }
        }
        if (num_splits) {
            /// We had to split the alignment; create a final subseg alignment
            aligns.push_back(x_CreateSubsegAlignment(last_seg, curr_seg - 1));
            ++num_splits;
        }
    }}
    break;

    case TSegs::e_Disc:
    {{
        const CSeq_align_set &segs = GetSegs().GetDisc();
        CSeq_align_set curr_disc;
        TSeqRange curr_range0;
        TSeqRange curr_range1;
        ITERATE (list< CRef<CSeq_align> >, seg_it, segs.Get()) {
            TSeqRange prev_range0 = curr_range0;
            TSeqRange prev_range1 = curr_range1;
            curr_range0 = (*seg_it)->GetSeqRange(0);
            curr_range1 = (*seg_it)->GetSeqRange(1);
            list< CRef<CSeq_align> > seg_splits;
            (*seg_it)->SplitOnLongDiscontinuity(seg_splits,
                                                discontinuity_threshold);

            if (prev_range0.Empty() ||
                (s_Distance(curr_range0, prev_range0) <= discontinuity_threshold &&
                 s_Distance(curr_range1, prev_range1) <= discontinuity_threshold))
            {
                /// distance between this alignment and previous one is less than
                /// threshold; first split section of this alignment should be joined
                /// to previous one
                curr_disc.Set().push_back(seg_splits.front());
                seg_splits.pop_front();
                if (seg_splits.empty()) {
                    continue;
                }
            }

            /// Put current contents of curr_disc in aligns
            aligns.push_back(s_GetJoinedAlignment(curr_disc));
            ++num_splits;

            /// Put results of current section's split in aligns except for last
            /// section; place last section in curr_disc, since it may have to be
            /// joined with next one
            num_splits += seg_splits.size();
            curr_disc.Set().clear();
            curr_disc.Set().push_back(seg_splits.back());
            seg_splits.pop_back();

            aligns.splice(aligns.end(), seg_splits);
        }
        if (num_splits) {
            /// We had to split the alignment; create a final alignment from
            /// the contents of curr_disc
            aligns.push_back(s_GetJoinedAlignment(curr_disc));
            ++num_splits;
        }
    }}
    break;

    default:
    break;
    }

    if (!num_splits) {
        aligns.push_back(CRef<CSeq_align>(const_cast<CSeq_align *>(this)));
    }
}

void CSeq_align::OffsetRow(TDim row,
                          TSignedSeqPos offset)
{
    if (offset == 0) return;

    switch (SetSegs().Which()) {
    case TSegs::e_Dendiag:
        NON_CONST_ITERATE(TSegs::TDendiag, dendiag_it, SetSegs().SetDendiag()) {
            (*dendiag_it)->OffsetRow(row, offset);
        }
        break;
    case TSegs::e_Denseg:
        SetSegs().SetDenseg().OffsetRow(row, offset);
        break;
    case TSegs::e_Std:
        NON_CONST_ITERATE(TSegs::TStd, std_it, SetSegs().SetStd()) {
            (*std_it)->OffsetRow(row, offset);
        }
        break;
    case TSegs::e_Disc:
        NON_CONST_ITERATE(CSeq_align_set::Tdata, seq_align_it, SetSegs().SetDisc().Set()) {
            (*seq_align_it)->OffsetRow(row, offset);
        }
        break;
    default:
        NCBI_THROW(CSeqalignException, eUnsupported,
                   "CSeq_align::OffsetRow() currently does not handle "
                   "this type of alignment");
    }
}


/// @deprecated
void CSeq_align::RemapToLoc(TDim row,
                            const CSeq_loc& dst_loc,
                            bool ignore_strand)
{
    // Limit to certain types of locs:
    switch (dst_loc.Which()) {
    case CSeq_loc::e_Whole:
        return;
    case CSeq_loc::e_Int:
        break;
    default:
        NCBI_THROW(CSeqalignException, eUnsupported,
                   "CSeq_align::RemapToLoc only supports int target seq-locs");
    }

    switch (SetSegs().Which()) {
    case TSegs::e_Denseg:
        SetSegs().SetDenseg().RemapToLoc(row, dst_loc, ignore_strand);
        break;
    case TSegs::e_Std:
        NON_CONST_ITERATE(TSegs::TStd, std_it, SetSegs().SetStd()) {
            (*std_it)->RemapToLoc(row, dst_loc, ignore_strand);
        }
        break;
    case TSegs::e_Disc:
        NON_CONST_ITERATE(CSeq_align_set::Tdata, seq_align_it, SetSegs().SetDisc().Set()) {
            (*seq_align_it)->RemapToLoc(row, dst_loc, ignore_strand);
        }
        break;
    default:
        NCBI_THROW(CSeqalignException, eUnsupported,
                   "CSeq_align::RemapToLoc only supports Dense-seg and Std-seg alignments.");
    }
}


CRef<CSeq_align> RemapAlignToLoc(const CSeq_align& align,
                                 CSeq_align::TDim  row,
                                 const CSeq_loc&   loc)
{
    if ( loc.IsWhole() ) {
        CRef<CSeq_align> copy(new CSeq_align);
        copy->Assign(align);
        return copy;
    }
    const CSeq_id* orig_id = loc.GetId();
    if ( !orig_id ) {
        NCBI_THROW(CAnnotMapperException, eBadLocation,
                   "Location with multiple ids can not be used to "
                   "remap seq-aligns.");
    }
    CRef<CSeq_id> id(new CSeq_id);
    id->Assign(*orig_id);

    // Create source seq-loc
    TSeqPos len = 0;
    for (CSeq_loc_CI it(loc); it; ++it) {
        if ( it.IsWhole() ) {
            NCBI_THROW(CAnnotMapperException, eBadLocation,
                    "Whole seq-loc can not be used to "
                    "remap seq-aligns.");
        }
        len += it.GetRange().GetLength();
    }
    CSeq_loc src_loc(*id, 0, len - 1);
    ENa_strand strand = loc.GetStrand();
    if (strand != eNa_strand_unknown) {
        src_loc.SetStrand(strand);
    }
    CSeq_loc_Mapper_Base mapper(src_loc, loc);
    return mapper.Map(align, row);
}

/// Get length of intersection between a range and a range collection
static inline TSeqPos s_IntersectionLength(const CRangeCollection<TSeqPos> &ranges,
                                           const TSeqRange &range)
{
    TSeqPos length = 0;
    ITERATE (CRangeCollection<TSeqPos>, it, ranges) {
        length += it->IntersectionWith(range).GetLength();
    }
    return length;
}

/// Retrieves the number of gaps in an alignment
/// @param align Alignment to examine [in]
/// @param get_total_count if true, the total number of gaps is retrived,
/// otherwise only the number of gap openings is retrieved
/// @throws CSeqalignException if alignment type is not supported
static TSeqPos
s_GetGapCount(const CSeq_align& align, CSeq_align::TDim row,
              const CRangeCollection<TSeqPos> &ranges,
              bool get_total_count)
{
    if (ranges.empty()) {
        return 0;
    }

    TSeqPos retval = 0;
    switch (align.GetSegs().Which()) {
    case CSeq_align::TSegs::e_Denseg:
        {{
            const CDense_seg& ds = align.GetSegs().GetDenseg();
            for (CDense_seg::TNumseg i = 0;  i < ds.GetNumseg();  ++i) {
                bool is_gapped = false;
                for (CDense_seg::TDim j = 0;  j < ds.GetDim();  ++j) {
                    if (ds.GetStarts()[i * ds.GetDim() + j] == -1 &&
                        (row < 0 || row == j))
                    {
                        is_gapped = true;
                        break;
                    }
                }
                if (is_gapped) {
                    TSeqPos gap_len = ds.GetLens()[i];
                    if (!ranges.begin()->IsWhole()) {
                        TSignedSeqPos gap_start = ds.GetStarts()[i * ds.GetDim()];
                        if (gap_start >= 0) {
                            gap_len = s_IntersectionLength(ranges,
                                TSeqRange(gap_start, gap_start + gap_len - 1));
                        } else {
                            gap_start = ds.GetStarts()[(i-1) * ds.GetDim()]
                                      + ds.GetLens()[i-1];
                            if (s_IntersectionLength(ranges,
                                    TSeqRange(gap_start, gap_start)) == 0)
                            {
                                gap_len = 0;
                            }
                        }
                    }
                    retval += (get_total_count ? gap_len : (gap_len ? 1 : 0));
                }
            }
        }}
        break;

    case CSeq_align::TSegs::e_Disc:
        {{
            ITERATE(CSeq_align::TSegs::TDisc::Tdata, iter,
                    align.GetSegs().GetDisc().Get()) {
                retval += s_GetGapCount(**iter, row, ranges, get_total_count);
            }
        }}
        break;

    case CSeq_align::TSegs::e_Spliced:
        {{
            ITERATE (CSpliced_seg::TExons, iter, align.GetSegs().GetSpliced().GetExons()) {
                for (CSeq_align::TDim r = 0; r < 2; ++r) {
                    if (row != r) {
                        CRangeCollection<TSeqPos> insertions =
                             (*iter)->GetRowSeq_insertions(
                                  r, align.GetSegs().GetSpliced(), ranges);
                        retval += get_total_count ? insertions.GetCoveredLength()
                                                  : insertions.size();
                    }
                }
            }
        }}
        break;

    case CSeq_align::TSegs::e_Std:
        if (row < 0 && !get_total_count && ranges.begin()->IsWhole()) {
            ITERATE (CSeq_align::TSegs::TStd, iter, align.GetSegs().GetStd()) {
                const CStd_seg& seg = **iter;
                ITERATE (CStd_seg::TLoc, it, seg.GetLoc()) {
                    const CSeq_loc& loc = **it;
                    if (loc.IsEmpty()) {
                        ++retval;
                        break;
                    }
                }
            }
            break;
        }

    default:
        NCBI_THROW(CSeqalignException, eUnsupported,
                   "CSeq_align::GetGapCount() currently does not handle "
                   "this type of alignment.");
    }

    return retval;
}

TSeqPos CSeq_align::GetTotalGapCount(TDim row) const
{
    return s_GetGapCount(*this, row,
                         CRangeCollection<TSeqPos>(TSeqRange::GetWhole()),
                         true);
}

TSeqPos CSeq_align::GetNumGapOpenings(TDim row) const
{
    return s_GetGapCount(*this, row,
                         CRangeCollection<TSeqPos>(TSeqRange::GetWhole()),
                         false);
}

TSeqPos CSeq_align::GetTotalGapCountWithinRange(
    const TSeqRange &range, TDim row) const
{
    return s_GetGapCount(*this, row, CRangeCollection<TSeqPos>(range), true);
}

TSeqPos CSeq_align::GetNumGapOpeningsWithinRange(
    const TSeqRange &range, TDim row) const
{
    return s_GetGapCount(*this, row, CRangeCollection<TSeqPos>(range), false);
}

TSeqPos CSeq_align::GetTotalGapCountWithinRanges(
    const CRangeCollection<TSeqPos> &ranges, TDim row) const
{
    return s_GetGapCount(*this, row, ranges, true);
}

TSeqPos CSeq_align::GetNumGapOpeningsWithinRanges(
    const CRangeCollection<TSeqPos> &ranges, TDim row) const
{
    return s_GetGapCount(*this, row, ranges, false);
}

string CSeq_align::SIndel::AsString() const
{
    return NStr::NumericToString(genomic_pos)
         + (row == 0 ? 'I' : 'D')
         + NStr::NumericToString(length);
}

static vector<CSeq_align::SIndel>
s_GetIndels(const CSeq_align& align, CSeq_align::TDim row,
            const CRangeCollection<TSeqPos> &ranges,
            bool include_frameshifts, bool include_non_frameshifts)
{
    vector<CSeq_align::SIndel> results;
    if (ranges.empty()) {
        return results;;
    }

    CRef<CSeq_align> generated_denseg;
    switch (align.GetSegs().Which()) {
    case CSeq_align::TSegs::e_Denseg:
        break;

    case CSeq_align::TSegs::e_Disc:
        generated_denseg = align.CreateDensegFromDisc();
        break;

    case CSeq_align::TSegs::e_Spliced:
        {{
            CRef<CSeq_align> as_disc_seg =
                align.GetSegs().GetSpliced().AsDiscSeg();
            if (align.GetSeqStart(1) > align.GetSeqStop(1)) {
                /// genomic sequence is circular; won't work on the entire
                /// alignment, need to check separately on the parts before and
                /// after the origin
                CSeq_align before_origin, after_origin;
                before_origin.SetType(CSeq_align::eType_disc);
                after_origin.SetType(CSeq_align::eType_disc);
                for (const CRef<CSeq_align> &segment
                     : as_disc_seg->GetSegs().GetDisc().Get())
                {
                    (segment->GetSeqStart(1) < align.GetSeqStart(1)
                       ? after_origin : before_origin)
                    . SetSegs().SetDisc().Set().push_back(segment);
                }
                vector<CSeq_align::SIndel> before_origin_frameshifts =
                    s_GetIndels(before_origin, row, ranges, include_frameshifts, include_non_frameshifts),
                                           after_origin_frameshifts =
                    s_GetIndels(after_origin, row, ranges, include_frameshifts, include_non_frameshifts);
                results = before_origin_frameshifts;
                results.insert(results.end(), after_origin_frameshifts.begin(),
                                              after_origin_frameshifts.end());
                return results;
            }
            generated_denseg = as_disc_seg->CreateDensegFromDisc();
        }}
        break;

    default:
        NCBI_THROW(CSeqalignException, eUnsupported,
                   "CSeq_align::GetNumFrameshifts() currently does not handle "
                   "this type of alignment.");
    }

    if (generated_denseg) {
        generated_denseg->SetSegs().SetDenseg().OrderAdjacentGaps();
        generated_denseg->SetSegs().SetDenseg().Compact();
        generated_denseg->SetSegs().SetDenseg().RemovePureGapSegs();
    }

    const CDense_seg &ds = (generated_denseg ? *generated_denseg : align)
                           . GetSegs().GetDenseg();

            for (CDense_seg::TNumseg i = 0;  i < ds.GetNumseg();  ++i) {
                bool is_gapped = false;
                for (CDense_seg::TDim j = 0;  j < ds.GetDim();  ++j) {
                    if (ds.GetStarts()[i * ds.GetDim() + j] == -1 && row != j)
                    {
                        is_gapped = true;
                        break;
                    }
                }
                if (is_gapped) {
                    TSeqPos gap_len = ds.GetLens()[i];
                    if (!ranges.begin()->IsWhole()) {
                        TSignedSeqPos gap_start = ds.GetStarts()[i * ds.GetDim()];
                        if (gap_start >= 0) {
                            gap_len = s_IntersectionLength(ranges,
                                TSeqRange(gap_start, gap_start + gap_len - 1));
                        } else {
                            gap_start = ds.GetStarts()[(i-1) * ds.GetDim()]
                                      + ds.GetLens()[i-1];
                            if (s_IntersectionLength(ranges,
                                    TSeqRange(gap_start, gap_start)) == 0)
                            {
                                gap_len = 0;
                            }
                        }
                    }
                    if (gap_len > 0 &&
                       ((include_frameshifts && gap_len % 3 != 0) ||
                        (include_non_frameshifts && gap_len % 3 == 0)))
                    {
                        TSignedSeqPos genomic_gap_start = ds.GetStarts()[i*ds.GetDim() + 1];
                        CSeq_align::TDim inserted_row = 1;
                        if (genomic_gap_start < 0) {
                            ENa_strand genomic_strand = ds.IsSetStrands() ? ds.GetStrands()[(i-1) * ds.GetDim() + 1] : eNa_strand_plus;
                            genomic_gap_start = genomic_strand == eNa_strand_minus
                                 ? ds.GetStarts()[(i+1) * ds.GetDim() + 1]
                                      + ds.GetLens()[i+1]
                                 : ds.GetStarts()[(i-1) * ds.GetDim() + 1]
                                      + ds.GetLens()[i-1];
                            inserted_row = 0;
                        }
                        results.push_back(CSeq_align::SIndel(genomic_gap_start, inserted_row, gap_len));
                    }
                }
            }

    return results;
}

TSeqPos CSeq_align::GetNumFrameshifts(TDim row) const
{
    return s_GetIndels(*this, row,
                         CRangeCollection<TSeqPos>(TSeqRange::GetWhole()), true, false)
           .size();
}

TSeqPos CSeq_align::GetNumFrameshiftsWithinRange(
    const TSeqRange &range, TDim row) const
{
    return s_GetIndels(*this, row, CRangeCollection<TSeqPos>(range), true, false)
           .size();
}

TSeqPos CSeq_align::GetNumFrameshiftsWithinRanges(
    const CRangeCollection<TSeqPos> &ranges, TDim row) const
{
    return s_GetIndels(*this, row, ranges, true, false)
           .size();
}

vector<CSeq_align::SIndel> CSeq_align::GetFrameshifts(TDim row) const
{
    return s_GetIndels(*this, row,
                         CRangeCollection<TSeqPos>(TSeqRange::GetWhole()), true, false);
}

vector<CSeq_align::SIndel> CSeq_align::GetFrameshiftsWithinRange(
    const TSeqRange &range, TDim row) const
{
    return s_GetIndels(*this, row, CRangeCollection<TSeqPos>(range), true, false);
}

vector<CSeq_align::SIndel> CSeq_align::GetFrameshiftsWithinRanges(
    const CRangeCollection<TSeqPos> &ranges, TDim row) const
{
    return s_GetIndels(*this, row, ranges, true, false);
}

vector<CSeq_align::SIndel> CSeq_align::GetNonFrameshifts(TDim row) const
{
    return s_GetIndels(*this, row,
                         CRangeCollection<TSeqPos>(TSeqRange::GetWhole()), false, true);
}

vector<CSeq_align::SIndel> CSeq_align::GetNonFrameshiftsWithinRange(
    const TSeqRange &range, TDim row) const
{
    return s_GetIndels(*this, row, CRangeCollection<TSeqPos>(range), false, true);
}

vector<CSeq_align::SIndel> CSeq_align::GetNonFrameshiftsWithinRanges(
    const CRangeCollection<TSeqPos> &ranges, TDim row) const
{
    return s_GetIndels(*this, row, ranges, false, true);
}

vector<CSeq_align::SIndel> CSeq_align::GetIndels(TDim row) const
{
    return s_GetIndels(*this, row,
                         CRangeCollection<TSeqPos>(TSeqRange::GetWhole()), true, true);
}

vector<CSeq_align::SIndel> CSeq_align::GetIndelsWithinRange(
    const TSeqRange &range, TDim row) const
{
    return s_GetIndels(*this, row, CRangeCollection<TSeqPos>(range), true, true);
}

vector<CSeq_align::SIndel> CSeq_align::GetIndelsWithinRanges(
    const CRangeCollection<TSeqPos> &ranges, TDim row) const
{
    return s_GetIndels(*this, row, ranges, true, true);
}


CRangeCollection<TSeqPos> CSeq_align::GetAlignedBases(TDim row) const
{
    CRangeCollection<TSeqPos> ranges;
    switch (GetSegs().Which()) {
    case CSeq_align::TSegs::e_Denseg:
        {{
            const CDense_seg& ds = GetSegs().GetDenseg();
            for (CDense_seg::TNumseg i = 0;  i < ds.GetNumseg();  ++i) {
                bool is_gapped = false;
                for (CDense_seg::TDim j = 0;  j < ds.GetDim();  ++j) {
                    if (ds.GetStarts()[i * ds.GetDim() + j] == -1)
                    {
                        is_gapped = true;
                        break;
                    }
                }
                if (!is_gapped) {
                    TSignedSeqPos start = ds.GetStarts()[i * ds.GetDim() + row];
                    ranges += TSeqRange(start, start + ds.GetLens()[i] - 1);
                }
            }
        }}
        break;

    case CSeq_align::TSegs::e_Disc:
        {{
            ITERATE(CSeq_align::TSegs::TDisc::Tdata, iter,
                    GetSegs().GetDisc().Get()) {
                ranges += (*iter)->GetAlignedBases(row);
            }
        }}
        break;

    case CSeq_align::TSegs::e_Spliced:
        {{
            ITERATE (CSpliced_seg::TExons, iter, GetSegs().GetSpliced().GetExons()) {
                ranges += (*iter)->GetRowSeq_range(row, true);
                ranges -= (*iter)->GetRowSeq_insertions(row, GetSegs().GetSpliced());
            }
        }}
        break;

    default:
        NCBI_THROW(CSeqalignException, eUnsupported,
                   "CSeq_align::GetInsertedRanges() currently does not handle "
                   "this type of alignment.");
    }

    return ranges;
}


/// Get length of dense-seg alignment within range
static inline TSeqPos s_DenseSegLength(const CDense_seg& ds,
                                       CDense_seg::TNumseg seg,
                                       const CRangeCollection<TSeqPos> &ranges)
{
    if(ranges.begin()->IsWhole())
        return ds.GetLens()[seg];

    TSignedSeqPos start = ds.GetStarts()[seg * ds.GetDim()];
    TSeqRange seg_range(start, start + ds.GetLens()[seg] - 1);
    return s_IntersectionLength(ranges, seg_range);
}

//////////////////////////////////////////////////////////////////////////////
///
/// calculate the length of our alignment within given range
///
static size_t s_GetAlignmentLength(const CSeq_align& align,
                                   const CRangeCollection<TSeqPos> &ranges,
                                   bool ungapped)
{
    if (ranges.empty()) {
        return 0;
    }

    size_t len = 0;
    switch (align.GetSegs().Which()) {
    case CSeq_align::TSegs::e_Denseg:
        {{
            const CDense_seg& ds = align.GetSegs().GetDenseg();
            if (ungapped) {
                for (CDense_seg::TNumseg i = 0;  i < ds.GetNumseg();  ++i) {
                    bool is_gap_seg = false;
                    for (CDense_seg::TDim j = 0;
                         !is_gap_seg  &&  j < ds.GetDim();  ++j) {
                        //int start = ds.GetStarts()[i * ds.GetDim() + j];
                        if (ds.GetStarts()[i * ds.GetDim() + j] == -1) {
                            is_gap_seg = true;
                        }
                    }
                    if ( !is_gap_seg ) {
                        len += s_DenseSegLength(ds, i, ranges);
                    }
                }
            } else {
                for (size_t i = 0;  i < ds.GetLens().size();  ++i) {
                    len += s_DenseSegLength(ds, i, ranges);
               }
            }
        }}
        break;

    case CSeq_align::TSegs::e_Disc:
        ITERATE (CSeq_align::TSegs::TDisc::Tdata, iter,
                 align.GetSegs().GetDisc().Get()) {
            len += s_GetAlignmentLength(**iter, ranges, ungapped);
        }
        break;

    case CSeq_align::TSegs::e_Std:
        {{
             if(!ranges.begin()->IsWhole())
                 NCBI_THROW(CSeqalignException, eUnsupported,
                            "Can't calculate alignment length within a range "
                            "for standard seg representation");

             /// pass 1:
             /// find total ranges
             vector<TSeqPos> sizes;

             size_t count = 0;
             ITERATE (CSeq_align::TSegs::TStd, iter, align.GetSegs().GetStd()) {
                 const CStd_seg& seg = **iter;

                 bool is_gap = false;
                 ITERATE (CStd_seg::TLoc, it, seg.GetLoc()) {
                     if ((*it)->IsEmpty()) {
                         is_gap = true;
                         break;
                     }
                 }

                 if ( !ungapped  ||  !is_gap) {
                     size_t i = 0;
                     ITERATE (CStd_seg::TLoc, it, seg.GetLoc()) {
                         if ( !(*it)->IsEmpty() ) {
                             if (sizes.empty()) {
                                 sizes.resize(seg.GetDim(), 0);
                             } else if (sizes.size() != (size_t)seg.GetDim()) {
                                 NCBI_THROW(CException, eUnknown,
                                            "invalid std-seg: "
                                            "inconsistent number of locs");
                             }

                             for (CSeq_loc_CI loc_it(**it);  loc_it;  ++loc_it) {
                                 sizes[i] += loc_it.GetRange().GetLength();
                             }
                         }

                         ++i;
                     }
                 }

                 ++count;
             }

             /// pass 2: determine shortest length
             if ( sizes.empty() ) return 0;
             vector<TSeqPos>::iterator iter = sizes.begin();
             vector<TSeqPos>::iterator smallest = iter;
             for (++iter;  iter != sizes.end();  ++iter) {
                 if (*iter < *smallest) {
                     smallest = iter;
                 }
             }
             return *smallest;
         }}
        break;

    case CSeq_align::TSegs::e_Spliced:
        ITERATE (CSpliced_seg::TExons, iter,
                 align.GetSegs().GetSpliced().GetExons()) {
            CRangeCollection<TSeqPos> exon_ranges = ranges;
            exon_ranges &= (*iter)->GetRowSeq_range(0, true);
            len += exon_ranges.GetCoveredLength();
            if (ungapped) {
                len -= (*iter)->GetRowSeq_insertions(
                         0, align.GetSegs().GetSpliced(), ranges).GetCoveredLength();
            } else {
                /// When computing length without gaps, genomic inserts
                /// are ignored; when computing with gaps, we need to add them
                /// to length if their starting point is within range
                len += (*iter)->GetRowSeq_insertions(
                         1, align.GetSegs().GetSpliced(), ranges).GetCoveredLength();
            }
        }
        break;

    default:
        _ASSERT(false);
        break;
    }

    return len;
}


TSeqPos CSeq_align::GetAlignLength(bool include_gaps) const
{
    return s_GetAlignmentLength(*this,
                                CRangeCollection<TSeqPos>(TSeqRange::GetWhole()),
                                !include_gaps );
}


TSeqPos CSeq_align::GetAlignLengthWithinRange(const TSeqRange &range,
                                              bool include_gaps) const
{
    return s_GetAlignmentLength(*this,
                                CRangeCollection<TSeqPos>(range),
                                !include_gaps );
}


TSeqPos CSeq_align::GetAlignLengthWithinRanges(const CRangeCollection<TSeqPos> &ranges,
                                               bool include_gaps) const
{
    return s_GetAlignmentLength(*this, ranges, !include_gaps );
}


double CSeq_align::AlignLengthRatio() const {
    TSeqRange r0 = GetSeqRange(0);
    TSeqRange r1 = GetSeqRange(1);
    double r = 0;
    if (r0.GetLength()) {
        r = double(r1.GetLength()) / double(r0.GetLength());
    }
    return r;
}


CRef<CSeq_loc> CSeq_align::CreateRowSeq_loc(TDim row) const
{
    CRef<CSeq_loc> loc(new CSeq_loc);
    switch (GetSegs().Which()) {
    case CSeq_align::TSegs::e_Dendiag:
        {
            ITERATE(TSegs::TDendiag, it, GetSegs().GetDendiag()) {
                loc->SetPacked_int().Set().push_back(
                    (*it)->CreateRowSeq_interval(row));
            }
            break;
        }
    case CSeq_align::TSegs::e_Denseg:
        {
            loc->SetInt(*GetSegs().GetDenseg().CreateRowSeq_interval(row));
            break;
        }
    case CSeq_align::TSegs::e_Std:
        {
            ITERATE(TSegs::TStd, it, GetSegs().GetStd()) {
                // Std-seg may contain empty locations, so
                // we have to use mix rather than packed-int.
                loc->SetMix().Set().push_back(
                    (*it)->CreateRowSeq_loc(row));
            }
            break;
        }
    case CSeq_align::TSegs::e_Disc:
        {
            ITERATE(TSegs::TDisc::Tdata, it, GetSegs().GetDisc().Get()) {
                loc->SetMix().Set().push_back((*it)->CreateRowSeq_loc(row));
            }
            break;
        }
    case CSeq_align::TSegs::e_Spliced:
        {
            if (row > 1) {
                NCBI_THROW(CSeqalignException, eInvalidRowNumber,
                           "CSeq_align::CreateRowSeq_loc() - row number must "
                           "be 0 or 1 for spliced-segs.");
            }
            const CSpliced_seg& spl = GetSegs().GetSpliced();
            ITERATE(TSegs::TSpliced::TExons, ex, spl.GetExons()) {
                // Spliced-seg may be required to get seq-id if
                // it's not set in the exon.
                loc->SetPacked_int().Set().push_back(
                    (*ex)->CreateRowSeq_interval(row, spl));
            }
            break;
        }
    case CSeq_align::TSegs::e_Packed:
    case CSeq_align::TSegs::e_Sparse:
    default:
        NCBI_THROW(CSeqalignException, eUnsupported,
                   "CSeq_align::CreateRowSeq_loc() currently does not handle "
                   "this type of alignment.");
        break;
    }
    return loc;
}

struct SLengthRange
{
    SLengthRange() : range(numeric_limits<TSeqPos>::max(), 0) {}
    void AddLength(TSeqPos length)
    {
        if (range.first > length) {
            range.first = length;
        }
        if (range.second < length) {
            range.second = length;
        }
    }
    void AddRange(CSeq_align::TLengthRange new_range)
    {
        if (range.first > new_range.first) {
            range.first = new_range.first;
        }
        if (range.second < new_range.second) {
            range.second = new_range.second;
        }
    }

    CSeq_align::TLengthRange range;
};

CSeq_align::TLengthRange CSeq_align::GapLengthRange() const
{
    SLengthRange length_range;
    switch (GetSegs().Which()) {
    case CSeq_align::TSegs::e_Denseg:
        {{
             const CDense_seg& ds = GetSegs().GetDenseg();
             for (CDense_seg::TNumseg i = 0;  i < ds.GetNumseg();  ++i) {
                 bool is_gapped = false;
                 for (CDense_seg::TDim j = 0;  j < ds.GetDim();  ++j) {
                     if (ds.GetStarts()[i * ds.GetDim() + j] == -1) {
                         is_gapped = true;
                         break;
                     }
                 }
                 if (is_gapped) {
                     length_range.AddLength(ds.GetLens()[i]);
                 }
             }
         }}
        break;

    case CSeq_align::TSegs::e_Disc:
        {{
             int num_rows = CheckNumRows();
             vector<TSeqRange> last_seg_ranges;
             ITERATE(CSeq_align::TSegs::TDisc::Tdata, iter,
                     GetSegs().GetDisc().Get()) {
                 length_range.AddRange((*iter)->GapLengthRange());
                 vector<TSeqRange> seg_ranges;
                 for (int i=0; i < num_rows; i++) {
                     seg_ranges.push_back((*iter)->GetSeqRange(i));
                     /// If this is not first segment, include gaps between last
                     /// segment and this one
                     if (!last_seg_ranges.empty()) {
                         TSeqPos gap = s_Distance(seg_ranges[i], last_seg_ranges[i]);
                         if (gap > 0) {
                             length_range.AddLength(gap);
                         }
                     }
                 }
                 last_seg_ranges = seg_ranges;
             }
         }}
        break;

    case CSeq_align::TSegs::e_Spliced:
        {{
             ITERATE (CSpliced_seg::TExons, iter,
                      GetSegs().GetSpliced().GetExons()) {
                 const CSpliced_exon& exon = **iter;
                 for (unsigned r = 0; r < 2; ++r) {
                     CRangeCollection<TSeqPos> insertions =
                          exon.GetRowSeq_insertions(
                               r, GetSegs().GetSpliced());
                     ITERATE (CRangeCollection<TSeqPos>, ins_it, insertions)
                     {
                         length_range.AddLength(ins_it->GetLength());
                     }
                 }
             }
         }}

        break;

    default:
        NCBI_THROW(CSeqalignException, eUnsupported,
                   "Can't get gap lengths for this type of alignment.");
    }
    return length_range.range;
}

CSeq_align::TLengthRange CSeq_align::IntronLengthRange() const
{
    if (!GetSegs().IsSpliced()) {
        NCBI_THROW(CSeqalignException, eUnsupported,
                   "Requested exon lengths for a non-spliced alignment.");
    }
    SLengthRange length_range;
    const CSpliced_exon* previous_exon = NULL;
    bool minus_strand = GetSeqStrand(1) == eNa_strand_minus;
    ITERATE (CSpliced_seg::TExons, iter,
             GetSegs().GetSpliced().GetExons()) {
        const CSpliced_exon* exon = &**iter;
        if (previous_exon) {
            bool descending_pos = exon->GetGenomic_end() <
                                previous_exon->GetGenomic_start();
            if (descending_pos != minus_strand) {
                /// Crossed origins; don't count as intron size
                continue;
            }
            TSeqRange intron((minus_strand ? exon : previous_exon)->GetGenomic_end() + 1,
                             (minus_strand ? previous_exon : exon)->GetGenomic_start() - 1);
            length_range.AddLength(intron.GetLength());
        }
        previous_exon = exon;
    }
    return length_range.range;
}

CSeq_align::TLengthRange CSeq_align::ExonLengthRange() const
{
    if (!GetSegs().IsSpliced()) {
        NCBI_THROW(CSeqalignException, eUnsupported,
                   "Requested exon lengths for a non-spliced alignment.");
    }
    SLengthRange length_range;
    ITERATE (CSpliced_seg::TExons, iter,
             GetSegs().GetSpliced().GetExons()) {
        const CSpliced_exon* exon = &**iter;
        length_range.AddLength(
            exon->GetGenomic_end() - exon->GetGenomic_start() + 1);
    }
    return length_range.range;
}

CRef<CSeq_align>  CSeq_align::x_CreateSubsegAlignment(int from, int to) const
{
    // break at this segment
    CRef<CSeq_align> align(new CSeq_align);
    align->SetDim(2);
    align->SetType(CSeq_align::eType_partial);
    const CDense_seg& seg = GetSegs().GetDenseg();
    CDense_seg& subseg = align->SetSegs().SetDenseg();

    subseg.SetIds() = seg.GetIds();
    subseg.SetDim(2);
    subseg.SetNumseg(to - from + 1);
    subseg.SetStarts().reserve(subseg.GetNumseg()*2);
    subseg.SetLens().reserve(subseg.GetNumseg());
    if (seg.IsSetStrands()) {
        subseg.SetStrands().reserve(subseg.GetNumseg()*2);
    }

    int i;

    for (i = from;  i <= to;  ++i) {
        subseg.SetLens().push_back(seg.GetLens()[i]);
        subseg.SetStarts().push_back(seg.GetStarts()[i * 2]);
        subseg.SetStarts().push_back(seg.GetStarts()[i * 2 + 1]);
        if (seg.IsSetStrands()) {
            subseg.SetStrands().push_back(seg.GetStrands()[i * 2]);
            subseg.SetStrands().push_back(seg.GetStrands()[i * 2 + 1]);
        }
    }

    subseg.TrimEndGaps();
    try {
        align->Validate(true);
    }
    catch (CException& e) {
        ERR_POST(Error << e);
        cerr << MSerial_AsnText << seg;
        cerr << MSerial_AsnText << *align;
        throw;
    }
    return align;
}


CConstRef<CUser_object> CSeq_align::FindExt(const string& ext_type) const
{
    CConstRef<CUser_object> ret;
    if ( IsSetExt() ) {
        ITERATE(TExt, it, GetExt()) {
            const CObject_id& obj_type = (*it)->GetType();
            if ( obj_type.IsStr()  &&  obj_type.GetStr() == ext_type ) {
                ret.Reset(it->GetPointer());
                break;
            }
        }
    }
    return ret;
}

CRef<CUser_object> CSeq_align::FindExt(const string& ext_type)
{
    CRef<CUser_object> ret;
    if ( IsSetExt() ) {
        NON_CONST_ITERATE(TExt, it, SetExt()) {
            const CObject_id& obj_type = (*it)->GetType();
            if ( obj_type.IsStr()  &&  obj_type.GetStr() == ext_type ) {
                ret.Reset(it->GetPointer());
                break;
            }
        }
    }
    return ret;
}

END_objects_SCOPE // namespace ncbi::objects::

END_NCBI_SCOPE