cav_alignset.cpp
上传用户:yhdzpy8989
上传日期:2007-06-13
资源大小:13604k
文件大小:15k
源码类别:

生物技术

开发平台:

C/C++

  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: cav_alignset.cpp,v $
  4.  * PRODUCTION Revision 1000.2  2004/06/01 19:41:12  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.5
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  9. /*  $Id: cav_alignset.cpp,v 1000.2 2004/06/01 19:41:12 gouriano Exp $
  10. * ===========================================================================
  11. *
  12. *                            PUBLIC DOMAIN NOTICE
  13. *               National Center for Biotechnology Information
  14. *
  15. *  This software/database is a "United States Government Work" under the
  16. *  terms of the United States Copyright Act.  It was written as part of
  17. *  the author's official duties as a United States Government employee and
  18. *  thus cannot be copyrighted.  This software/database is freely available
  19. *  to the public for use. The National Library of Medicine and the U.S.
  20. *  Government have not placed any restriction on its use or reproduction.
  21. *
  22. *  Although all reasonable efforts have been taken to ensure the accuracy
  23. *  and reliability of the software and data, the NLM and the U.S.
  24. *  Government do not and cannot warrant the performance or results that
  25. *  may be obtained by using this software or data. The NLM and the U.S.
  26. *  Government disclaim all warranties, express or implied, including
  27. *  warranties of performance, merchantability or fitness for any particular
  28. *  purpose.
  29. *
  30. *  Please cite the author in any work or product based on this material.
  31. *
  32. * ===========================================================================
  33. *
  34. * Authors:  Paul Thiessen
  35. *
  36. * File Description:
  37. *      Classes to hold sets of master/slave pairwise alignments
  38. *
  39. * ===========================================================================
  40. */
  41. #include <ncbi_pch.hpp>
  42. #include <objects/seqalign/Seq_align.hpp>
  43. #include <objects/seqalign/Dense_diag.hpp>
  44. #include <objects/seqalign/Dense_seg.hpp>
  45. #include <objects/seqloc/Seq_id.hpp>
  46. #include <objects/seqloc/PDB_seq_id.hpp>
  47. #include <objects/general/Object_id.hpp>
  48. #include <objects/seqloc/Textseq_id.hpp>
  49. #include <objects/seqloc/PDB_mol_id.hpp>
  50. #include <objects/seq/Bioseq.hpp>
  51. #include <map>
  52. #include <objtools/cddalignview/cav_alignset.hpp>
  53. #include <objtools/cddalignview/cav_seqset.hpp>
  54. #include <objtools/cddalignview/cddalignview.h>
  55. BEGIN_NCBI_SCOPE
  56. USING_SCOPE(objects);
  57. typedef list < const CSeq_align * > SeqAlignList;
  58. typedef vector < CRef < CSeq_id > > SeqIdList;
  59. static bool IsAMatch(const Sequence *seq, const CSeq_id& sid)
  60. {
  61.     if (sid.IsGi()) {
  62.         if (sid.GetGi() == seq->gi)
  63.             return true;
  64.         return false;
  65.     }
  66.     if (sid.IsPdb()) {
  67.         if (sid.GetPdb().GetMol().Get() == seq->pdbID) {
  68.             if (sid.GetPdb().IsSetChain()) {
  69.                 if (sid.GetPdb().GetChain() != seq->pdbChain) {
  70.                     return false;
  71.                 }
  72.             }
  73.             return true;
  74.         }
  75.         return false;
  76.     }
  77.     if (sid.IsLocal() && sid.GetLocal().IsStr()) {
  78.         if (sid.GetLocal().GetStr() == seq->pdbID)
  79.             return true;
  80.         return false;
  81.     }
  82.     if (sid.IsGenbank() && sid.GetGenbank().IsSetAccession()) {
  83.         if (sid.GetGenbank().GetAccession() == seq->accession)
  84.             return true;
  85.         return false;
  86.     }
  87.     if (sid.IsSwissprot() && sid.GetSwissprot().IsSetAccession()) {
  88.         return (sid.GetSwissprot().GetAccession() == seq->accession);
  89.     }
  90.     ERR_POST(Error << "IsAMatch - can't match this type of Seq-id");
  91.     return false;
  92. }
  93. AlignmentSet::AlignmentSet(SequenceSet *sequenceSet, const SeqAnnotList& seqAnnots) :
  94.     master(NULL), status(CAV_ERROR_ALIGNMENTS)
  95. {
  96.     SeqAlignList seqaligns;
  97.     // first, make a list of alignments
  98.     SeqAnnotList::const_iterator n, ne = seqAnnots.end();
  99.     for (n=seqAnnots.begin(); n!=ne; ++n) {
  100.         if (!n->GetObject().GetData().IsAlign()) {
  101.             ERR_POST(Error << "AlignmentSet::AlignmentSet() - confused by alignment data format");
  102.             return;
  103.         }
  104.         CSeq_annot::C_Data::TAlign::const_iterator
  105.             a, ae = n->GetObject().GetData().GetAlign().end();
  106.         for (a=n->GetObject().GetData().GetAlign().begin(); a!=ae; ++a) {
  107.             // verify this is a type of alignment we can deal with
  108.             if (!(a->GetObject().GetType() != CSeq_align::eType_partial ||
  109.                   a->GetObject().GetType() != CSeq_align::eType_diags) ||
  110.                 !a->GetObject().IsSetDim() || a->GetObject().GetDim() != 2 ||
  111.                 (!a->GetObject().GetSegs().IsDendiag() && !a->GetObject().GetSegs().IsDenseg())) {
  112.                 ERR_POST(Error << "AlignmentSet::AlignmentSet() - confused by alignment type");
  113.                 return;
  114.             }
  115.             seqaligns.push_back(a->GetPointer());
  116.         }
  117.     }
  118.     if (seqaligns.size() == 0) {
  119.         ERR_POST(Error << "AlignmentSet::AlignmentSet() - no valid Seq-aligns present");
  120.         return;
  121.     }
  122.     // figure out which sequence is the master
  123.     // if there aren't any structures, then the first sequence of each pair must be the same,
  124.     // and is the master
  125.     else {
  126.         SeqAlignList::const_iterator a, ae = seqaligns.end();
  127.         for (a=seqaligns.begin(); a!=ae; ++a) {
  128.             const SeqIdList& sids = (*a)->GetSegs().IsDendiag() ?
  129.                 (*a)->GetSegs().GetDendiag().front()->GetIds() :
  130.                 (*a)->GetSegs().GetDenseg().GetIds();
  131.             if (!master) {
  132.                 SequenceSet::SequenceList::const_iterator
  133.                     s, se = sequenceSet->sequences.end();
  134.                 for (s=sequenceSet->sequences.begin(); s!=se; ++s) {
  135.                     if (IsAMatch(*s, sids.front().GetObject())) {
  136.                         master = *s;
  137.                         break;
  138.                     }
  139.                 }
  140.             }
  141.             if (!IsAMatch(master, sids.front().GetObject())) {
  142.                 ERR_POST(Error << "AlignmentSet::AlignmentSet() - master must be first sequence of every alignment");
  143.                 return;
  144.             }
  145.         }
  146.         sequenceSet->master = master;
  147.     }
  148.     if (master) {
  149.         ERR_POST(Info << "determined that master sequence is " << master->GetTitle());
  150.     } else {
  151.         ERR_POST(Error << "AlignmentSet::AlignmentSet() - couldn't determine which is master sequence");
  152.         return;
  153.     }
  154.     // then make an alignment from each Seq-align; for any sequence that has structure, make
  155.     // sure that that Sequence object only appears once in the alignment
  156.     SeqAlignList::const_iterator s, se = seqaligns.end();
  157.     int i = 1;
  158.     for (s=seqaligns.begin(); s!=se; ++s, ++i) {
  159.         const MasterSlaveAlignment *alignment =
  160.             new MasterSlaveAlignment(sequenceSet, master, **s);
  161.         if (!alignment || alignment->Status() != CAV_SUCCESS) {
  162.             ERR_POST(Error << "Error parsing master/slave pairwise alignment #" << i);
  163.             status = alignment->Status();
  164.             return;
  165.         }
  166.         alignments.push_back(alignment);
  167.     }
  168.     ERR_POST(Info << "number of alignments: " << alignments.size());
  169.     status = CAV_SUCCESS;
  170. }
  171. AlignmentSet::~AlignmentSet(void)
  172. {
  173.     for (int i=0; i<alignments.size(); ++i) delete alignments[i];
  174. }
  175. ///// MasterSlaveAlignment methods /////
  176. static bool SameSequence(const Sequence *s, const CBioseq& bioseq)
  177. {
  178.     CBioseq::TId::const_iterator i, ie = bioseq.GetId().end();
  179.     for (i=bioseq.GetId().begin(); i!=ie; ++i) {
  180.         if (IsAMatch(s, **i))
  181.             return true;
  182.     }
  183.     return false;
  184. }
  185. MasterSlaveAlignment::MasterSlaveAlignment(
  186.     const SequenceSet *sequenceSet,
  187.     const Sequence *masterSequence,
  188.     const CSeq_align& seqAlign) :
  189.     master(masterSequence), slave(NULL), status(CAV_ERROR_PAIRWISE)
  190. {
  191.     // resize alignment and block vector
  192.     masterToSlave.resize(master->Length(), -1);
  193.     SequenceSet::SequenceList::const_iterator
  194.         s = sequenceSet->sequences.begin(),
  195.         se = sequenceSet->sequences.end();
  196.     // find slave sequence for this alignment, and order (master or slave first)
  197.     const SeqIdList& sids = seqAlign.GetSegs().IsDendiag() ?
  198.         seqAlign.GetSegs().GetDendiag().front()->GetIds() :
  199.         seqAlign.GetSegs().GetDenseg().GetIds();
  200.     bool masterFirst = true;
  201.     for (; s!=se; ++s) {
  202.         if (SameSequence(*s, master->bioseqASN.GetObject())) continue;
  203.         if (IsAMatch(*s, sids.back().GetObject())) {
  204.             break;
  205.         } else if (IsAMatch(*s, sids.front().GetObject())) {
  206.             masterFirst = false;
  207.             break;
  208.         }
  209.     }
  210.     if (s == se) {
  211.         // special case of master seq. aligned to itself
  212.         if (IsAMatch(master, sids.back().GetObject()) && IsAMatch(master, sids.front().GetObject())) {
  213.             slave = master;
  214.         } else {
  215.             ERR_POST(Error << "MasterSlaveAlignment::MasterSlaveAlignment() - n"
  216.                 "couldn't find matching unaligned slave sequence");
  217.             return;
  218.         }
  219.     } else {
  220.         slave = *s;
  221.     }
  222.     int i, masterRes, slaveRes;
  223.     // unpack dendiag alignment
  224.     if (seqAlign.GetSegs().IsDendiag()) {
  225.         CSeq_align::C_Segs::TDendiag::const_iterator d , de = seqAlign.GetSegs().GetDendiag().end();
  226.         for (d=seqAlign.GetSegs().GetDendiag().begin(); d!=de; ++d) {
  227.             const CDense_diag& block = d->GetObject();
  228.             if (!block.IsSetDim() || block.GetDim() != 2 ||
  229.                 block.GetIds().size() != 2 ||
  230.                 block.GetStarts().size() != 2) {
  231.                 ERR_POST(Error << "MasterSlaveAlignment::MasterSlaveAlignment() - n"
  232.                     "incorrect dendiag block dimensions");
  233.                 return;
  234.             }
  235.             // make sure identities of master and slave sequences match in each block
  236.             if ((masterFirst &&
  237.                     (!IsAMatch(master, block.GetIds().front().GetObject()) ||
  238.                      !IsAMatch(slave, block.GetIds().back().GetObject()))) ||
  239.                 (!masterFirst &&
  240.                     (!IsAMatch(master, block.GetIds().back().GetObject()) ||
  241.                      !IsAMatch(slave, block.GetIds().front().GetObject())))) {
  242.                 ERR_POST(Error << "MasterSlaveAlignment::MasterSlaveAlignment() - n"
  243.                     "mismatched Seq-id in dendiag block");
  244.                 return;
  245.             }
  246.             // finally, actually unpack the data into the alignment vector
  247.             for (i=0; i<block.GetLen(); ++i) {
  248.                 if (masterFirst) {
  249.                     masterRes = block.GetStarts().front() + i;
  250.                     slaveRes = block.GetStarts().back() + i;
  251.                 } else {
  252.                     masterRes = block.GetStarts().back() + i;
  253.                     slaveRes = block.GetStarts().front() + i;
  254.                 }
  255.                 if (masterRes >= master->Length() || slaveRes >= slave->Length()) {
  256.                     ERR_POST(Critical << "MasterSlaveAlignment::MasterSlaveAlignment() - n"
  257.                         "seqloc in dendiag block > length of sequence!");
  258.                     return;
  259.                 }
  260.                 masterToSlave[masterRes] = slaveRes;
  261.             }
  262.         }
  263.     }
  264.     // unpack denseg alignment
  265.     else if (seqAlign.GetSegs().IsDenseg()) {
  266.         const CDense_seg& block = seqAlign.GetSegs().GetDenseg();
  267.         if (!block.IsSetDim() && block.GetDim() != 2 ||
  268.             block.GetIds().size() != 2 ||
  269.             block.GetStarts().size() != 2 * block.GetNumseg() ||
  270.             block.GetLens().size() != block.GetNumseg()) {
  271.             ERR_POST(Error << "MasterSlaveAlignment::MasterSlaveAlignment() - n"
  272.                 "incorrect denseg block dimension");
  273.             return;
  274.         }
  275.         // make sure identities of master and slave sequences match in each block
  276.         if ((masterFirst &&
  277.                 (!IsAMatch(master, block.GetIds().front().GetObject()) ||
  278.                  !IsAMatch(slave, block.GetIds().back().GetObject()))) ||
  279.             (!masterFirst &&
  280.                 (!IsAMatch(master, block.GetIds().back().GetObject()) ||
  281.                  !IsAMatch(slave, block.GetIds().front().GetObject())))) {
  282.             ERR_POST(Error << "MasterSlaveAlignment::MasterSlaveAlignment() - n"
  283.                 "mismatched Seq-id in denseg block");
  284.             return;
  285.         }
  286.         // finally, actually unpack the data into the alignment vector
  287.         CDense_seg::TStarts::const_iterator starts = block.GetStarts().begin();
  288.         CDense_seg::TLens::const_iterator lens, le = block.GetLens().end();
  289.         for (lens=block.GetLens().begin(); lens!=le; ++lens) {
  290.             if (masterFirst) {
  291.                 masterRes = *(starts++);
  292.                 slaveRes = *(starts++);
  293.             } else {
  294.                 slaveRes = *(starts++);
  295.                 masterRes = *(starts++);
  296.             }
  297.             if (masterRes != -1 && slaveRes != -1) { // skip gaps
  298.                 if ((masterRes + *lens - 1) >= master->Length() ||
  299.                     (slaveRes + *lens - 1) >= slave->Length()) {
  300.                     ERR_POST(Critical << "MasterSlaveAlignment::MasterSlaveAlignment() - n"
  301.                         "seqloc in denseg block > length of sequence!");
  302.                     return;
  303.                 }
  304.                 for (i=0; i<*lens; ++i)
  305.                     masterToSlave[masterRes + i] = slaveRes + i;
  306.             }
  307.         }
  308.     }
  309.     status = CAV_SUCCESS;
  310. }
  311. END_NCBI_SCOPE
  312. /*
  313. * ---------------------------------------------------------------------------
  314. * $Log: cav_alignset.cpp,v $
  315. * Revision 1000.2  2004/06/01 19:41:12  gouriano
  316. * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.5
  317. *
  318. * Revision 1.5  2004/05/21 21:42:51  gorelenk
  319. * Added PCH ncbi_pch.hpp
  320. *
  321. * Revision 1.4  2004/05/08 10:56:20  thiessen
  322. * better handling of repeated sequences
  323. *
  324. * Revision 1.3  2004/03/15 18:51:27  thiessen
  325. * prefer prefix vs. postfix ++/-- operators
  326. *
  327. * Revision 1.2  2003/06/02 16:06:41  dicuccio
  328. * Rearranged src/objects/ subtree.  This includes the following shifts:
  329. *     - src/objects/asn2asn --> arc/app/asn2asn
  330. *     - src/objects/testmedline --> src/objects/ncbimime/test
  331. *     - src/objects/objmgr --> src/objmgr
  332. *     - src/objects/util --> src/objmgr/util
  333. *     - src/objects/alnmgr --> src/objtools/alnmgr
  334. *     - src/objects/flat --> src/objtools/flat
  335. *     - src/objects/validator --> src/objtools/validator
  336. *     - src/objects/cddalignview --> src/objtools/cddalignview
  337. * In addition, libseq now includes six of the objects/seq... libs, and libmmdb
  338. * replaces the three libmmdb? libs.
  339. *
  340. * Revision 1.1  2003/03/19 19:04:12  thiessen
  341. * move again
  342. *
  343. * Revision 1.1  2003/03/19 05:33:43  thiessen
  344. * move to src/app/cddalignview
  345. *
  346. * Revision 1.7  2003/02/03 17:52:03  thiessen
  347. * move CVS Log to end of file
  348. *
  349. * Revision 1.6  2002/07/11 17:16:54  thiessen
  350. * change for STL type in object lib
  351. *
  352. * Revision 1.5  2002/02/08 19:53:17  thiessen
  353. * add annotation to text/HTML displays
  354. *
  355. * Revision 1.4  2001/07/12 21:32:12  thiessen
  356. * update for swissprot accessions
  357. *
  358. * Revision 1.3  2001/01/29 18:13:33  thiessen
  359. * split into C-callable library + main
  360. *
  361. * Revision 1.2  2001/01/22 15:55:11  thiessen
  362. * correctly set up ncbi namespacing
  363. *
  364. * Revision 1.1  2001/01/22 13:15:23  thiessen
  365. * initial checkin
  366. *
  367. */