开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 数值算法/人工智能 > 生物技术 > 查看源码
validerror_align.cpp
资源名称：ncbi_cxx_win--Jul_08_2004.zip　[点击查看]
上传用户：yhdzpy8989
上传日期：2007-06-13
资源大小：13604k
文件大小：35k
源码类别：

生物技术

开发平台：

C/C++

validerror_align.cpp：源码内容
  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: validerror_align.cpp,v $
  4.  * PRODUCTION Revision 1000.1  2004/06/01 19:47:45  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.8
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  10. /*  $Id: validerror_align.cpp,v 1000.1 2004/06/01 19:47:45 gouriano Exp $
  11.  * ===========================================================================
  12.  *
  13.  *                            PUBLIC DOMAIN NOTICE
  14.  *               National Center for Biotechnology Information
  15.  *
  16.  *  This software/database is a "United States Government Work" under the
  17.  *  terms of the United States Copyright Act.  It was written as part of
  18.  *  the author's official duties as a United States Government employee and
  19.  *  thus cannot be copyrighted.  This software/database is freely available
  20.  *  to the public for use. The National Library of Medicine and the U.S.
  21.  *  Government have not placed any restriction on its use or reproduction.
  22.  *
  23.  *  Although all reasonable efforts have been taken to ensure the accuracy
  24.  *  and reliability of the software and data, the NLM and the U.S.
  25.  *  Government do not and cannot warrant the performance or results that
  26.  *  may be obtained by using this software or data. The NLM and the U.S.
  27.  *  Government disclaim all warranties, express or implied, including
  28.  *  warranties of performance, merchantability or fitness for any particular
  29.  *  purpose.
  30.  *
  31.  *  Please cite the author in any work or product based on this material.
  32.  *
  33.  * ===========================================================================
  34.  *
  35.  * Author:  Jonathan Kans, Clifford Clausen, Aaron Ucko......
  36.  *
  37.  * File Description:
  38.  *   validation of seq_align
  39.  *   .......
  40.  *
  41.  */
  42. #include <ncbi_pch.hpp>
  43. #include <corelib/ncbistd.hpp>
  45. #include <objects/seqalign/Seq_align.hpp>
  46. #include <objects/seqalign/Dense_diag.hpp>
  47. #include <objects/seqalign/Dense_seg.hpp>
  48. #include <objects/seqalign/Packed_seg.hpp>
  49. #include <objects/seqalign/Seq_align_set.hpp>
  50. #include <objects/seqalign/Std_seg.hpp>
  51. #include <objmgr/util/sequence.hpp>
  53. #include <map>
  54. #include <vector>
  55. #include <algorithm>
  57. #include "validatorp.hpp"
  60. BEGIN_NCBI_SCOPE
  61. BEGIN_SCOPE(objects)
  62. BEGIN_SCOPE(validator)
  64. USING_SCOPE(sequence);
  67. // ================================  Public  ================================
  70. CValidError_align::CValidError_align(CValidError_imp& imp) :
  71.     CValidError_base(imp)
  72. {
  73. }
  76. CValidError_align::~CValidError_align(void)
  77. {
  78. }
  81. void CValidError_align::ValidateSeqAlign(const CSeq_align& align)
  82. {
  83.     const CSeq_align::TSegs& segs = align.GetSegs();
  84.     CSeq_align::C_Segs::E_Choice segtype = segs.Which();
  85.     switch ( segtype ) {
  86.     
  87.     case CSeq_align::C_Segs::e_Dendiag:
  88.         x_ValidateDendiag(segs.GetDendiag(), align);
  89.         break;
  91.     case CSeq_align::C_Segs::e_Denseg:
  92.         x_ValidateDenseg(segs.GetDenseg(), align);
  93.         break;
  95.     case CSeq_align::C_Segs::e_Std:
  96.         x_ValidateStd(segs.GetStd(), align);
  97.         break;
  98.     case CSeq_align::C_Segs::e_Packed:
  99.         x_ValidatePacked(segs.GetPacked(), align);
  100.         break;
  102.     case CSeq_align::C_Segs::e_Disc:
  103.         // call recursively
  104.         ITERATE(CSeq_align_set::Tdata, sali, segs.GetDisc().Get()) {
  105.             ValidateSeqAlign(**sali);
  106.         }
  107.         return;
  109.     case CSeq_align::C_Segs::e_not_set:
  110.     default:
  111.         PostErr(eDiag_Error, eErr_SEQ_ALIGN_Segtype,
  112.             "This alignment has undefined or unsupported Seq_align"
  113.             "segtype " + NStr::IntToString(segtype) + ".", align);
  114.         break;
  115.     }  // end of switch statement
  116. }
  119. // ================================  Private  ===============================
  122. void CValidError_align::x_ValidateDenseg
  123. (const TDenseg& denseg,
  124.  const CSeq_align& align)
  125. {
  126.     bool cont = true;
  128.     // assert dim >= 2
  129.     if ( !x_ValidateDim(denseg, align) ) {
  130.         return;
  131.     }
  133.     size_t dim     = denseg.GetDim();
  134.     size_t numseg  = denseg.GetNumseg();
  136.     // assert dim == Ids.size()
  137.     if ( dim != denseg.GetIds().size() ) {
  138.         PostErr(eDiag_Error, eErr_SEQ_ALIGN_SegsDimMismatch,
  139.                 "Mismatch between specified dimension (" +
  140.                 NStr::UIntToString(dim) + 
  141.                 ") and number of Seq-ids (" + 
  142.                 NStr::UIntToString(denseg.GetIds().size()) + ")",
  143.                 align);
  144.         cont = false;
  145.     }
  147.     // assert numseg == Lens.size()
  148.     if ( numseg != denseg.GetLens().size() ) {
  149.         PostErr(eDiag_Error, eErr_SEQ_ALIGN_SegsNumsegMismatch,
  150.                 "Mismatch between specified numseg (" +
  151.                 NStr::UIntToString(numseg) + 
  152.                 ") and number of Lens (" + 
  153.                 NStr::UIntToString(denseg.GetLens().size()) + ")",
  154.                 align);
  155.         cont = false;
  156.     }
  158.     // assert dim * numseg == Starts.size()
  159.     if ( dim * numseg != denseg.GetStarts().size() ) {
  160.         PostErr(eDiag_Error, eErr_SEQ_ALIGN_SegsStartsMismatch,
  161.                 "The number of Starts (" + 
  162.                 NStr::UIntToString(denseg.GetStarts().size()) +
  163.                 ") does not match the expected size of dim * numseg (" +
  164.                 NStr::UIntToString(dim * numseg) + ")", align);
  165.         cont = false;
  166.     } 
  168.     if( cont ) {
  169.         x_ValidateStrand(denseg, align);
  170.         x_ValidateFastaLike(denseg, align);
  171.         x_ValidateSegmentGap(denseg, align);
  172.         
  173.         if ( m_Imp.IsRemoteFetch() ) {
  174.             x_ValidateSeqId(align);
  175.             x_ValidateSeqLength(denseg, align);
  176.         }
  177.     }
  178. }
  183. void CValidError_align::x_ValidatePacked
  184. (const TPacked& packed,
  185.  const CSeq_align& align)
  186. {
  187.     bool cont = true;
  189.     // assert dim >= 2
  190.     if ( !x_ValidateDim(packed, align) ) {
  191.         return;
  192.     }
  194.     size_t dim     = packed.GetDim();
  195.     size_t numseg  = packed.GetNumseg();
  196.     
  197.     // assert dim == Ids.size()
  198.     if ( dim != packed.GetIds().size() ) {
  199.         PostErr(eDiag_Error, eErr_SEQ_ALIGN_SegsDimMismatch,
  200.             "Mismatch between specified dimension (" +
  201.             NStr::UIntToString(dim) + 
  202.             ") and number of Seq-ids (" + 
  203.             NStr::UIntToString(packed.GetIds().size()) + ")",
  204.             align);
  205.         cont = false;
  206.     }
  207.     
  208.     // assert numseg == Lens.size()
  209.     if ( numseg != packed.GetLens().size() ) {
  210.         PostErr(eDiag_Error, eErr_SEQ_ALIGN_SegsDimMismatch,
  211.             "Mismatch between specified numseg (" +
  212.             NStr::UIntToString(numseg) + 
  213.             ") and number of Lens (" + 
  214.             NStr::UIntToString(packed.GetLens().size()) + ")",
  215.             align);
  216.         cont = false;
  217.     }
  218.     
  219.     // assert dim * numseg == Present.size()
  220.     if ( dim * numseg != packed.GetPresent().size() ) {
  221.         PostErr(eDiag_Error, eErr_SEQ_ALIGN_SegsPresentMismatch,
  222.                 "The number of Present (" + 
  223.                 NStr::UIntToString(packed.GetPresent().size()) +
  224.                 ") does not match the expected size of dim * numseg (" +
  225.                 NStr::UIntToString(dim * numseg) + ")", align);
  226.         cont = false;
  227.     } else {
  228.         // assert # of '1' bits in present == Starts.size()
  229.         size_t bits = x_CountBits(packed.GetPresent());
  230.         if ( bits != packed.GetStarts().size() ) {
  231.             PostErr(eDiag_Error, eErr_SEQ_ALIGN_SegsPresentStartsMismatch,
  232.                 "Starts does not have the same number of elements (" +
  233.                 NStr::UIntToString(packed.GetStarts().size()) + 
  234.                 ") as specified by Present (" + NStr::UIntToString(bits) +
  235.                 ")", align);
  236.             cont = false;
  237.         }
  239.         // assert # of '1' bits in present == Strands.size() (if exists)
  240.         if ( packed.IsSetStrands() ) {
  241.             if ( bits != packed.GetStarts().size() ) {
  242.                 PostErr(eDiag_Error, eErr_SEQ_ALIGN_SegsPresentStrandsMismatch,
  243.                     "Strands does not have the same number of elements (" +
  244.                     NStr::UIntToString(packed.GetStrands().size()) +
  245.                     ") as specified by Present (" +
  246.                     NStr::UIntToString(bits) + ")", align);
  247.                 cont = false;
  248.             }
  249.         }
  250.     }
  251.     
  252.     if( cont ) {
  253.         x_ValidateStrand(packed, align);
  254.         x_ValidateFastaLike(packed, align);
  255.         x_ValidateSegmentGap(packed, align);
  256.         
  257.         if ( m_Imp.IsRemoteFetch() ) {
  258.             x_ValidateSeqId(align);
  259.             x_ValidateSeqLength(packed, align);
  260.         }
  261.     }
  262. }
  265. size_t CValidError_align::x_CountBits(const CPacked_seg::TPresent& present)
  266. {
  267.     size_t bits = 0;
  268.     ITERATE( CPacked_seg::TPresent, iter, present ) {
  269.         if ( *iter != 0 ) {
  270.             Uchar mask = 0x01;
  271.             for ( size_t i = 0; i < 7; ++i ) {
  272.                 bits += (*iter & mask) ? 1 : 0;
  273.                 mask = mask << 1;
  274.             }
  275.         }
  276.     }
  277.     return bits;
  278. }
  281. void CValidError_align::x_ValidateDendiag
  282. (const TDendiag& dendiags,
  283.  const CSeq_align& align)
  284. {
  285.     size_t num_dendiag = 0;
  286.     ITERATE( TDendiag, dendiag_iter, dendiags ) {
  287.         bool cont = true;
  288.         ++num_dendiag;
  290.         const CDense_diag& dendiag = **dendiag_iter;
  291.         size_t dim = dendiag.GetDim();
  293.         // assert dim >= 2
  294.         if ( !x_ValidateDim(dendiag, align, num_dendiag) ) {
  295.             continue;
  296.         }
  298.         // assert dim == Ids.size()
  299.         if ( dim != dendiag.GetIds().size() ) {
  300.             PostErr(eDiag_Error, eErr_SEQ_ALIGN_SegsDimMismatch,
  301.                 "Mismatch between specified dimension (" +
  302.                 NStr::UIntToString(dim) + 
  303.                 ") and number of Seq-ids (" + 
  304.                 NStr::UIntToString(dendiag.GetIds().size()) + ") in dendiag " +
  305.                 NStr::UIntToString(num_dendiag), align);
  306.             cont = false;
  307.         }
  309.         // assert dim == Starts.size()
  310.         if ( dim != dendiag.GetStarts().size() ) {
  311.             PostErr(eDiag_Error, eErr_SEQ_ALIGN_SegsDimMismatch,
  312.                 "Mismatch between specified dimension (" +
  313.                 NStr::UIntToString(dim) + 
  314.                 ") and number ofStarts (" + 
  315.                 NStr::UIntToString(dendiag.GetStarts().size()) + 
  316.                 ") in dendiag " + NStr::UIntToString(num_dendiag), align);
  317.             cont = false;
  318.         } 
  320.         // assert dim == Strands.size() (if exist)
  321.         if ( dendiag.IsSetStrands() ) {
  322.             if ( dim != dendiag.GetStrands().size() ) {
  323.                 PostErr(eDiag_Error, eErr_SEQ_ALIGN_SegsDimMismatch,
  324.                     "Mismatch between specified dimension (" +
  325.                     NStr::UIntToString(dim) + 
  326.                     ") and number of Strands (" + 
  327.                     NStr::UIntToString(dendiag.GetStrands().size()) + 
  328.                     ") in dendiag " + NStr::UIntToString(num_dendiag), align);
  329.                 cont = false;
  330.             } 
  331.         }
  333.         if ( cont ) {
  334.             if ( m_Imp.IsRemoteFetch() ) {
  335.                 x_ValidateSeqId(align);
  336.                 x_ValidateSeqLength(dendiag, num_dendiag, align);
  337.             }
  338.         }
  339.     }
  340. }
  343. void CValidError_align::x_ValidateStd
  344. (const TStd& std_segs,
  345.  const CSeq_align& align)
  346. {
  347.     bool cont = true;
  349.     size_t num_stdseg = 0;
  350.     ITERATE( TStd, stdseg_iter, std_segs) {
  351.         ++num_stdseg;
  353.         const CStd_seg& stdseg = **stdseg_iter;
  354.         size_t dim = stdseg.GetDim();
  356.         // assert dim >= 2
  357.         if ( !x_ValidateDim(stdseg, align, num_stdseg) ) {
  358.             cont = false;
  359.             continue;
  360.         }
  362.         // assert dim == Loc.size()
  363.         if ( dim != stdseg.GetLoc().size() ) {
  364.             PostErr(eDiag_Error, eErr_SEQ_ALIGN_SegsDimMismatch,
  365.                 "Mismatch between specified dimension (" +
  366.                 NStr::UIntToString(dim) + 
  367.                 ") and number of Seq-locs (" + 
  368.                 NStr::UIntToString(stdseg.GetLoc().size()) + ") in stdseg " +
  369.                 NStr::UIntToString(num_stdseg), align);
  370.             cont = false;
  371.             continue;
  372.         }
  374.         // assert dim == Ids.size()
  375.         if ( stdseg.IsSetIds() ) {
  376.             if ( dim != stdseg.GetIds().size() ) {
  377.                 PostErr(eDiag_Error, eErr_SEQ_ALIGN_SegsDimMismatch,
  378.                     "Mismatch between specified dimension (" +
  379.                     NStr::UIntToString(dim) + 
  380.                     ") and number of Seq-ids (" + 
  381.                     NStr::UIntToString(stdseg.GetIds().size()) + ")",
  382.                     align);
  383.                 cont = false;
  384.                 continue;
  385.             }
  386.         }
  387.     }
  389.     if( cont ) {
  390.         x_ValidateStrand(std_segs, align);
  391.         x_ValidateFastaLike(std_segs, align);
  392.         x_ValidateSegmentGap(std_segs, align);
  393.         
  394.         if ( m_Imp.IsRemoteFetch() ) {
  395.             x_ValidateSeqId(align);
  396.             x_ValidateSeqLength(std_segs, align);
  397.         }
  398.     }
  399. }
  402. template <typename T>
  403. bool CValidError_align::x_ValidateDim
  404. (T& obj,
  405.  const CSeq_align& align,
  406.  size_t part)
  407. {
  408.     if ( !obj.IsSetDim() ) {
  409.         string msg = "Dim not set";
  410.         if ( part > 0 ) {
  411.             msg += " Part " + NStr::UIntToString(part);
  412.         }
  413.         PostErr(eDiag_Error, eErr_SEQ_ALIGN_SegsInvalidDim, msg, align);
  414.         return false;
  415.     }
  417.     size_t dim = obj.GetDim();        
  418.     if ( dim < 2 ) {
  419.         string msg = "Dimension (" + NStr::UIntToString(dim) + 
  420.             ") must be at least 2.";
  421.         if ( part > 0 ) {
  422.             msg += " Part " + NStr::UIntToString(part);
  423.         }
  424.         PostErr(eDiag_Error, eErr_SEQ_ALIGN_SegsInvalidDim, msg, align);
  425.         return false;
  426.     }
  427.     return true;
  428. }
  431. //===========================================================================
  432. // x_ValidateStrand:
  433. //
  434. //  Check if the  strand is consistent in SeqAlignment of global
  435. //  or partial type.
  436. //===========================================================================
  438. void CValidError_align::x_ValidateStrand
  439. (const TDenseg& denseg,
  440.  const CSeq_align& align)
  441. {
  442.     if ( !denseg.IsSetStrands() ) {
  443.         return;
  444.     }
  445.     
  446.     size_t dim = denseg.GetDim();
  447.     size_t numseg = denseg.GetNumseg();
  448.     const CDense_seg::TStrands& strands = denseg.GetStrands();
  450.     // go through id for each alignment sequence
  451.     for ( size_t id = 0; id < dim; ++id ) {
  452.         ENa_strand strand1 = strands[id];
  453.         
  454.         for ( size_t seg = 0; seg < numseg; ++seg ) {
  455.             ENa_strand strand2 = strands[id + (seg * dim)];
  457.             // skip undefined strand
  458.             if ( strand2 == eNa_strand_unknown  ||
  459.                  strand2 == eNa_strand_other ) {
  460.                 continue;
  461.             }
  463.             if ( strand1 == eNa_strand_unknown  ||
  464.                  strand1 == eNa_strand_other ) {
  465.                 strand1 = strand2;
  466.                 continue;
  467.             }
  469.             // strands should be same for a given seq-id
  470.             if ( strand1 != strand2 ) {
  471.                 PostErr(eDiag_Error, eErr_SEQ_ALIGN_StrandRev,
  472.                     "The strand labels for SeqId " + 
  473.                     denseg.GetIds()[id]->AsFastaString() +
  474.                     " are inconsistent across the alignment. "
  475.                     "The first inconsistent region is the " + 
  476.                     NStr::UIntToString(seg), align);
  477.                     break;
  478.             }
  479.         }
  480.     }
  481. }
  484. void CValidError_align::x_ValidateStrand
  485. (const TPacked& packed,
  486.  const CSeq_align& align)
  487. {
  488.     if ( !packed.IsSetStrands() ) {
  489.         return;
  490.     }
  492.     size_t dim = packed.GetDim();
  493.     const CPacked_seg::TPresent& present = packed.GetPresent();
  494.     const CPacked_seg::TStrands& strands = packed.GetStrands();
  495.     CPacked_seg::TStrands::const_iterator strand = strands.begin();
  496.     
  497.     vector<ENa_strand> id_strands(dim, eNa_strand_unknown);
  499.     Uchar mask = 0x01;
  500.     size_t present_size = present.size();
  501.     for ( size_t i = 0; i < present_size; ++i ) {
  502.         if ( present[i] == 0 ) {   // no bits in this char
  503.             continue;
  504.         }
  505.         for ( int shift = 7; shift >= 0; --shift ) {
  506.             if ( present[i] & (mask << shift) ) {
  507.                 size_t id = i * 8 + 7 - shift;
  508.                 if ( *strand == eNa_strand_unknown  ||  
  509.                      *strand == eNa_strand_other ) {
  510.                     ++strand;
  511.                     continue;
  512.                 }
  513.                 
  514.                 if ( id_strands[id] == eNa_strand_unknown  || 
  515.                      id_strands[id] == eNa_strand_other ) {
  516.                     id_strands[id] = *strand;
  517.                     ++strand;
  518.                     continue;
  519.                 }
  520.                 
  521.                 if ( id_strands[id] != *strand ) {
  522.                     CPacked_seg::TIds::const_iterator id_iter =
  523.                         packed.GetIds().begin();
  524.                     for ( ; id; --id ) {
  525.                         ++id_iter;
  526.                     }
  527.                     
  528.                     PostErr(eDiag_Error, eErr_SEQ_ALIGN_StrandRev,
  529.                         "Strand for SeqId " + (*id_iter)->AsFastaString() +
  530.                         " are inconsistent across the alignment", align);
  531.                 }
  532.                 ++strand;
  533.             }
  534.         }
  535.     }
  536. }
  539. void CValidError_align::x_ValidateStrand
  540. (const TStd& std_segs,
  541.  const CSeq_align& align)
  542. {
  543.     map< CConstRef<CSeq_id>, ENa_strand > strands;
  545.     ITERATE ( TStd, stdseg, std_segs ) {
  546.         ITERATE ( CStd_seg::TLoc, loc_iter, (*stdseg)->GetLoc() ) {
  547.             const CSeq_loc& loc = **loc_iter;
  548.             
  549.             if ( !IsOneBioseq(loc, m_Scope) ) {
  550.                 // !!! should probably be an error
  551.                 continue;
  552.             }
  553.             CConstRef<CSeq_id> id(&GetId(loc, m_Scope));
  554.             ENa_strand strand = GetStrand(loc, m_Scope);
  556.             if ( strand == eNa_strand_unknown  || 
  557.                  strand == eNa_strand_other ) {
  558.                 continue;
  559.             }
  561.             if ( strands[id] == eNa_strand_unknown  ||
  562.                  strands[id] == eNa_strand_other ) {
  563.                 strands[id] = strand;
  564.             }
  566.             if ( strands[id] != strand ) {
  567.                 PostErr(eDiag_Error, eErr_SEQ_ALIGN_StrandRev,
  568.                     "Strands for SeqId " + id->AsFastaString() + 
  569.                     " are inconsistent across the alignment", align);
  570.             }
  571.         }
  572.     }
  573. }
  576. //===========================================================================
  577. // x_ValidateFastaLike:
  578. //
  579. //  Check if an alignment is FASTA-like. 
  580. //  Alignment is FASTA-like if all gaps are at the end with dimensions > 2.
  581. //===========================================================================
  583. void CValidError_align::x_ValidateFastaLike
  584. (const TDenseg& denseg,
  585.  const CSeq_align& align)
  586. {
  587.     // check only global or partial type
  588.     if ( (align.GetType() != CSeq_align::eType_global  &&
  589.          align.GetType() != CSeq_align::eType_partial) ||
  590.          denseg.GetDim() <= 2 ) {
  591.         return;
  592.     }
  594.     size_t dim = denseg.GetDim();
  595.     size_t numseg = denseg.GetNumseg();
  597.     vector<string> fasta_like;
  599.     for ( size_t id = 0; id < dim; ++id ) {
  600.         bool gap = false;
  601.         
  602.         const CDense_seg::TStarts& starts = denseg.GetStarts();
  603.         for ( size_t seg = 0; seg < numseg; ++ seg ) {
  604.             // if start value is -1, set gap flag to true
  605.             if ( starts[id + (dim * seg)] < 0 ) {
  606.                 gap = true;
  607.             } else if ( gap  ) {
  608.                 // if a positive start value is found after the initial -1 
  609.                 // start value, it's not fasta like. 
  610.                 //no need to check this sequence further
  611.                 break;
  612.             } 
  614.             if ( seg == numseg - 1  &&  gap ) {
  615.                 // if no more positive start value are found after the initial
  616.                 // -1 start value, it's fasta like
  617.                 fasta_like.push_back(denseg.GetIds()[id]->AsFastaString());
  618.             }
  619.         }
  620.     }
  622.     if ( !fasta_like.empty() ) {
  623.         string fasta_like_ids;
  624.         bool first = true;
  625.         ITERATE( vector<string>, idstr, fasta_like ) {
  626.             if ( first ) {
  627.                 first = false;
  628.             } else {
  629.                 fasta_like_ids += ", ";
  630.             }
  632.             fasta_like_ids += *idstr;
  633.         }
  634.         fasta_like_ids += ".";
  636.         PostErr(eDiag_Error, eErr_SEQ_ALIGN_FastaLike,
  637.             "This may be a fasta-like alignment for SeqIds: " + 
  638.             fasta_like_ids, align);
  639.     }
  641. }           
  644. void CValidError_align::x_ValidateFastaLike
  645. (const TPacked& packed,
  646.  const CSeq_align& align)
  647. {
  648.     // check only global or partial type
  649.     if ( align.GetType() == CSeq_align::eType_global  ||
  650.          align.GetType() == CSeq_align::eType_partial ||
  651.          packed.GetDim() <= 2 ) {
  652.         return;
  653.     }
  655.     static Uchar bits[] = { 0x80, 0x40, 0x20, 0x10, 0x08, 0x04, 0x02, 0x01 };
  656.     
  657.     size_t dim = packed.GetDim();
  658.     size_t numseg = packed.GetNumseg();
  660.     CPacked_seg::TIds::const_iterator id_iter = packed.GetIds().begin();
  661.     const CPacked_seg::TPresent& present = packed.GetPresent();
  663.     for ( size_t id = 0; id < dim;  ++id, ++id_iter ) {
  664.         bool gap = false;
  665.         for ( size_t seg = 0; seg < numseg; ++ seg ) {
  666.             // if a segment is not present, set gap flag to true
  667.             size_t i = id + (dim * seg);
  668.             if ( !(present[i / 8] & bits[i % 8]) ) {
  669.                 gap = true;
  670.             } else if ( gap ) {
  671.                 // if a segment is found later, then it's not fasta like, 
  672.                 // no need to check this sequence further.
  673.                 break;
  674.             } 
  676.             if ( seg == numseg - 1  &&  gap ) {
  677.                 // if no more segments are found for this sequence, 
  678.                 // then it's  fasta like.
  679.                 PostErr(eDiag_Error, eErr_SEQ_ALIGN_FastaLike,
  680.                     "This may be a fasta-like alignment for SeqId: " + 
  681.                     (*id_iter)->AsFastaString(), align);
  682.             }
  683.         }
  684.     }
  685. }
  688. void CValidError_align::x_ValidateFastaLike
  689. (const TStd& std_segs,
  690.  const CSeq_align& align)
  691. {
  692.     // check only global or partial type
  693.     if ( align.GetType() == CSeq_align::eType_global  ||
  694.          align.GetType() == CSeq_align::eType_partial ) {
  695.         return;
  696.     }
  697.     // suspected seq-ids
  698.     typedef map< CConstRef<CSeq_id>, bool > TIdGapMap;
  699.     TIdGapMap fasta_like;
  701.     ITERATE ( TStd, stdseg, std_segs ) {
  702.         ITERATE ( CStd_seg::TLoc, loc_iter, (*stdseg)->GetLoc() ) {
  703.             const CSeq_loc& loc = **loc_iter;
  704.             
  705.             if ( !IsOneBioseq(loc, m_Scope) ) {
  706.                 // !!! should probably be an error
  707.                 continue;
  708.             }
  709.             CConstRef<CSeq_id> id(&GetId(loc, m_Scope));
  710.             bool gap = loc.IsEmpty()  ||  loc.IsNull();
  711.             if ( gap  &&  fasta_like.find(id) == fasta_like.end() ) {
  712.                 fasta_like[id] = true;
  713.             }
  714.             if ( !gap  &&  fasta_like.find(id) != fasta_like.end() ) {
  715.                 fasta_like[id] = false;
  716.             }
  717.         }
  718.     }
  719.     
  720.     ITERATE ( TIdGapMap, iter, fasta_like ) {
  721.         if ( iter->second ) {
  722.             PostErr(eDiag_Error, eErr_SEQ_ALIGN_FastaLike,
  723.                 "This may be a fasta-like alignment for SeqId: " + 
  724.                 iter->first->AsFastaString(), align);
  725.         }
  726.     }
  727. }
  730. //===========================================================================
  731. // x_ValidateSegmentGap:
  732. //
  733. // Check if there is a gap for all sequences in a segment.
  734. //===========================================================================
  736. void CValidError_align::x_ValidateSegmentGap
  737. (const TDenseg& denseg,
  738.  const CSeq_align& align)
  739. {
  740.     int numseg  = denseg.GetNumseg();
  741.     int dim     = denseg.GetDim();
  742.     const CDense_seg::TStarts& starts = denseg.GetStarts();
  744.     for ( int seg = 0; seg < numseg; ++seg ) {
  745.         bool seggap = true;
  746.         for ( int id = 0; id < dim; ++id ) {
  747.             if ( starts[seg * dim + id] != -1 ) {
  748.                 seggap = false;
  749.                 break;
  750.             }
  751.         }
  752.         if ( seggap ) {
  753.             // no sequence is present in this segment
  754.             PostErr(eDiag_Error, eErr_SEQ_ALIGN_SegmentGap,
  755.                 "Segment " + NStr::UIntToString(seg) + " contains only gaps.",
  756.                 align);
  757.         }
  758.     }
  759. }
  762. void CValidError_align::x_ValidateSegmentGap
  763. (const TPacked& packed,
  764.  const CSeq_align& align)
  765. {
  766.     static Uchar bits[] = { 0x80, 0x40, 0x20, 0x10, 0x08, 0x04, 0x02, 0x01 };
  768.     size_t numseg = packed.GetNumseg();
  769.     size_t dim = packed.GetDim();
  770.     const CPacked_seg::TPresent& present = packed.GetPresent();
  772.     for ( size_t seg = 0; seg < numseg;  ++seg) {
  773.         size_t id = 0;
  774.         for ( ; id < dim; ++id ) {
  775.             size_t i = id + (dim * seg);
  776.             if ( (present[i / 8] & bits[i % 8]) ) {
  777.                 break;
  778.             }
  779.         }
  780.         if ( id == dim ) {
  781.             // no sequence is present in this segment
  782.             PostErr(eDiag_Error, eErr_SEQ_ALIGN_SegmentGap,
  783.                 "Segment " + NStr::UIntToString(seg) + "contains only gaps.",
  784.                 align);
  785.         }
  786.     }       
  787. }
  790. void CValidError_align::x_ValidateSegmentGap
  791. (const TStd& std_segs,
  792.  const CSeq_align& align)
  793. {
  794.     size_t seg = 0;
  795.     ITERATE ( TStd, stdseg, std_segs ) {
  796.         bool gap = true;
  797.         ITERATE ( CStd_seg::TLoc, loc, (*stdseg)->GetLoc() ) {
  798.             if ( !(*loc)->IsEmpty()  ||  !(*loc)->IsEmpty() ) {
  799.                 gap = false;
  800.                 break;
  801.             }
  802.         }
  803.         if ( gap ) {
  804.             // no sequence is present in this segment
  805.             PostErr(eDiag_Error, eErr_SEQ_ALIGN_SegmentGap,
  806.                 "Segment " + NStr::UIntToString(seg) + "contains only gaps.",
  807.                 align);
  808.         }
  809.         ++seg;
  810.     }
  811. }
  814. //===========================================================================
  815. // x_ValidateSeqIdInSeqAlign:
  816. //
  817. //  Validate SeqId in sequence alignment.
  818. //===========================================================================
  820. void CValidError_align::x_ValidateSeqId(const CSeq_align& align)
  821. {
  822.     vector< CRef< CSeq_id > > ids;
  823.     x_GetIds(align, ids);
  825.     ITERATE( vector< CRef< CSeq_id > >, id_iter, ids ) {
  826.         const CSeq_id& id = **id_iter;
  827.         if ( id.IsLocal() ) {
  828.             if ( !m_Scope->GetBioseqHandle(id) ) {
  829.                 PostErr(eDiag_Error, eErr_SEQ_ALIGN_SeqIdProblem,
  830.                     "The sequence corresponding to SeqId " + 
  831.                     id.AsFastaString() + " could not be found.",
  832.                     align);
  833.             }
  834.         }
  835.     }
  836. }
  839. void CValidError_align::x_GetIds
  840. (const CSeq_align& align,
  841.  vector< CRef< CSeq_id > >& ids)
  842. {
  843.     ids.clear();
  845.     switch ( align.GetSegs().Which() ) {
  847.     case CSeq_align::C_Segs::e_Dendiag:
  848.         ITERATE( TDendiag, diag_seg, align.GetSegs().GetDendiag() ) {
  849.             const vector< CRef< CSeq_id > >& diag_ids = (*diag_seg)->GetIds();
  850.             copy(diag_ids.begin(), diag_ids.end(), back_inserter(ids));
  851.         }
  852.         break;
  853.         
  854.     case CSeq_align::C_Segs::e_Denseg:
  855.         ids = align.GetSegs().GetDenseg().GetIds();
  856.         break;
  857.         
  858.     case CSeq_align::C_Segs::e_Packed:
  859.         copy(align.GetSegs().GetPacked().GetIds().begin(),
  860.              align.GetSegs().GetPacked().GetIds().end(),
  861.              back_inserter(ids));
  862.         break;
  863.         
  864.     case CSeq_align::C_Segs::e_Std:
  865.         ITERATE( TStd, std_seg, align.GetSegs().GetStd() ) {
  866.             ITERATE( CStd_seg::TLoc, loc, (*std_seg)->GetLoc() ) {
  867.                 CSeq_id* idp = const_cast<CSeq_id*>(&GetId(**loc, m_Scope));
  868.                 CRef<CSeq_id> ref(idp);
  869.                 ids.push_back(ref);
  870.             }
  871.         }
  872.         break;
  873.             
  874.     default:
  875.         break;
  876.     }
  877. }
  880. //===========================================================================
  881. // x_ValidateSeqLength:
  882. //
  883. //  Check segment length, start and end point in Dense_diag, Dense_seg,  
  884. //  Packed_seg and Std_seg.
  885. //===========================================================================
  887. // Make sure that, in Dense_diag alignment, segment length is not greater
  888. // than Bioseq length
  889. void CValidError_align::x_ValidateSeqLength
  890. (const CDense_diag& dendiag,
  891.  size_t dendiag_num,
  892.  const CSeq_align& align)
  893. {
  894.     size_t dim = dendiag.GetDim();
  895.     TSeqPos len = dendiag.GetLen();
  896.     const CDense_diag::TIds& ids = dendiag.GetIds();
  897.     
  898.     CDense_diag::TStarts::const_iterator starts_iter = 
  899.             dendiag.GetStarts().begin();
  900.     
  901.     for ( size_t id = 0; id < dim; ++id ) {
  902.         TSeqPos bslen = GetLength(*(ids[id]), m_Scope);
  903.         TSeqPos start = *starts_iter;
  905.         // verify start
  906.         if ( start > bslen ) {
  907.             PostErr(eDiag_Error, eErr_SEQ_ALIGN_StartMorethanBiolen,
  908.                     "Start (" + NStr::UIntToString(start) +
  909.                     ") exceeds bioseq length (" +
  910.                     NStr::UIntToString(bslen) +
  911.                     ") for seq-id " + ids[id]->AsFastaString() +
  912.                     "in dendiag " + NStr::UIntToString(dendiag_num),
  913.                     align);
  914.         }
  915.         
  916.         // verify length
  917.         if ( start + len > bslen ) {
  918.             PostErr(eDiag_Error, eErr_SEQ_ALIGN_SumLenStart,
  919.                     "Start + length (" + NStr::UIntToString(start + len) +
  920.                     ") exceeds bioseq length (" +
  921.                     NStr::UIntToString(bslen) +
  922.                     ") for seq-id " + ids[id]->AsFastaString() +
  923.                     "in dendiag " + NStr::UIntToString(dendiag_num),
  924.                     align);
  925.         }
  926.         ++starts_iter;
  927.     }
  928. }
  931.         
  933. void CValidError_align::x_ValidateSeqLength
  934. (const TDenseg& denseg,
  935.  const CSeq_align& align)
  936. {
  937.     int dim     = denseg.GetDim();
  938.     int numseg  = denseg.GetNumseg();
  939.     const CDense_seg::TIds& ids       = denseg.GetIds();
  940.     const CDense_seg::TStarts& starts = denseg.GetStarts();
  941.     const CDense_seg::TLens& lens      = denseg.GetLens();
  942.     bool minus = false;
  944.     for ( int id = 0; id < dim; ++id ) {
  945.         TSeqPos bslen = GetLength(*(ids[id]), m_Scope);
  946.         minus = denseg.IsSetStrands()  &&
  947.             denseg.GetStrands()[id] == eNa_strand_minus;
  948.         
  949.         for ( int seg = 0; seg < numseg; ++seg ) {
  950.             size_t curr_index = 
  951.                 id + (minus ? numseg - seg - 1 : seg) * dim;
  952.             // no need to verify if segment is not present
  953.             if ( starts[curr_index] == -1 ) {
  954.                 continue;
  955.             }
  956.             size_t lens_index = minus ? numseg - seg - 1 : seg;
  958.             // verify that start plus segment does not exceed total bioseq len
  959.             if ( starts[curr_index] + lens[lens_index] > bslen ) {
  960.                 PostErr(eDiag_Error, eErr_SEQ_ALIGN_SumLenStart,
  961.                     "Start + segment length (" + 
  962.                     NStr::UIntToString(starts[curr_index] + lens[lens_index]) +
  963.                     ") exceeds bioseq length (" +
  964.                     NStr::UIntToString(bslen) + ")", align);
  965.             }
  967.             // find the next segment that is present
  968.             size_t next_index = curr_index;
  969.             int next_seg;
  970.             for ( next_seg = seg + 1; next_seg < numseg; ++next_seg ) {
  971.                 next_index = 
  972.                     id + (minus ? numseg - next_seg - 1 : next_seg) * dim;
  973.                 
  974.                 if ( starts[next_index] != -1 ) {
  975.                     break;
  976.                 }
  977.             }
  978.             if ( next_seg == numseg  ||  next_index == curr_index ) {
  979.                 continue;
  980.             }
  982.             // length plus start should be equal to the closest next 
  983.             // start that is not -1
  984.             if ( starts[curr_index] + (TSignedSeqPos)lens[lens_index] !=
  985.                 starts[next_index] ) {
  986.                 PostErr(eDiag_Error, eErr_SEQ_ALIGN_DensegLenStart,
  987.                     "Start + segment length (" + 
  988.                     NStr::UIntToString(starts[curr_index] + lens[lens_index]) +
  989.                     ") exceeds next start (" +
  990.                     NStr::UIntToString(starts[next_index]) + ")", align);
  991.             }
  992.         }
  993.     }
  994. }
  997. void CValidError_align::x_ValidateSeqLength
  998. (const TPacked& packed,
  999.  const CSeq_align& align)
  1000. {
  1001.     static Uchar bits[] = { 0x80, 0x40, 0x20, 0x10, 0x08, 0x04, 0x02, 0x01 };
  1003.     size_t dim = packed.GetDim();
  1004.     size_t numseg = packed.GetNumseg();
  1005.     //CPacked_seg::TStarts::const_iterator start =
  1006.     //    packed.GetStarts().begin();
  1008.     const CPacked_seg::TPresent& present = packed.GetPresent();
  1010.     for ( size_t id = 0; id < dim; ++id ) {
  1011.         for ( size_t seg = 0; seg < numseg; ++seg ) {
  1012.             size_t i = id + seg * dim;
  1013.             if ( (present[i / 8] & bits[i % 8]) ) {
  1014.                 // !!!
  1015.             }
  1016.         }
  1017.     }
  1018. }
  1021. void CValidError_align::x_ValidateSeqLength
  1022. (const TStd& std_segs,
  1023.  const CSeq_align& align)
  1024. {
  1025.     ITERATE( TStd, iter, std_segs ) {
  1026.         const CStd_seg& stdseg = **iter;
  1028.         ITERATE ( CStd_seg::TLoc, loc_iter, stdseg.GetLoc() ) {
  1029.             const CSeq_loc& loc = **loc_iter;
  1030.     
  1031.             if ( loc.IsWhole()  || loc.IsEmpty()  ||  loc.IsNull() ) {
  1032.                 continue;
  1033.             }
  1035.             if ( !IsOneBioseq(loc, m_Scope) ) {
  1036.                 continue;
  1037.             }
  1039.             TSeqPos from = loc.GetTotalRange().GetFrom();
  1040.             TSeqPos to   = loc.GetTotalRange().GetTo();
  1041.             TSeqPos loclen = GetLength( loc, m_Scope);
  1042.             TSeqPos bslen = GetLength(GetId(loc, m_Scope), m_Scope);
  1043.             string  bslen_str = NStr::UIntToString(bslen);
  1044.             string label;
  1045.             loc.GetLabel(&label);
  1047.             if ( from > bslen - 1 ) { 
  1048.                 PostErr(eDiag_Error, eErr_SEQ_ALIGN_StartMorethanBiolen,
  1049.                     "Loaction: " + label + ". From (" + 
  1050.                     NStr::UIntToString(from) + 
  1051.                     ") is more than bioseq length (" + bslen_str + ")", 
  1052.                     align);
  1053.             }
  1055.             if ( to > bslen - 1 ) { 
  1056.                 PostErr(eDiag_Error, eErr_SEQ_ALIGN_EndMorethanBiolen,
  1057.                     "Loaction: " + label + ". To (" + NStr::UIntToString(to) +
  1058.                     ") is more than bioseq length (" + bslen_str + ")", align);
  1059.             }
  1061.             if ( loclen > bslen ) {
  1062.                 PostErr(eDiag_Error, eErr_SEQ_ALIGN_LenMorethanBiolen,
  1063.                     "Loaction: " + label + ". Length (" + 
  1064.                     NStr::UIntToString(loclen) + 
  1065.                     ") is more than bioseq length (" + bslen_str + ")", align);
  1066.             }
  1067.         }
  1068.     }
  1069. }
  1072. END_SCOPE(validator)
  1073. END_SCOPE(objects)
  1074. END_NCBI_SCOPE
  1077. /*
  1078. * ===========================================================================
  1079. *
  1080. * $Log: validerror_align.cpp,v $
  1081. * Revision 1000.1  2004/06/01 19:47:45  gouriano
  1082. * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.8
  1083. *
  1084. * Revision 1.8  2004/05/21 21:42:56  gorelenk
  1085. * Added PCH ncbi_pch.hpp
  1086. *
  1087. * Revision 1.7  2003/09/30 18:30:38  shomrat
  1088. * changed unsigned to signed to eliminate infinite loop
  1089. *
  1090. * Revision 1.6  2003/06/02 16:06:43  dicuccio
  1091. * Rearranged src/objects/ subtree.  This includes the following shifts:
  1092. *     - src/objects/asn2asn --> arc/app/asn2asn
  1093. *     - src/objects/testmedline --> src/objects/ncbimime/test
  1094. *     - src/objects/objmgr --> src/objmgr
  1095. *     - src/objects/util --> src/objmgr/util
  1096. *     - src/objects/alnmgr --> src/objtools/alnmgr
  1097. *     - src/objects/flat --> src/objtools/flat
  1098. *     - src/objects/validator --> src/objtools/validator
  1099. *     - src/objects/cddalignview --> src/objtools/cddalignview
  1100. * In addition, libseq now includes six of the objects/seq... libs, and libmmdb
  1101. * replaces the three libmmdb? libs.
  1102. *
  1103. * Revision 1.5  2003/05/28 16:24:40  shomrat
  1104. * Report a single FastaLike error from each alignment
  1105. *
  1106. * Revision 1.4  2003/04/29 14:58:07  shomrat
  1107. * Implemented SeqAlign validation
  1108. *
  1109. * Revision 1.3  2003/03/31 14:40:52  shomrat
  1110. * $id: -> $id$
  1111. *
  1112. * Revision 1.2  2002/12/24 16:52:42  shomrat
  1113. * Changes to include directives
  1114. *
  1115. * Revision 1.1  2002/12/23 20:15:59  shomrat
  1116. * Initial submission after splitting former implementation
  1117. *
  1118. *
  1119. * ===========================================================================
  1120. */