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

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

生物技术

开发平台：

C/C++

blast_seqalign.cpp：源码内容
  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: blast_seqalign.cpp,v $
  4.  * PRODUCTION Revision 1000.6  2004/06/01 18:05:52  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.42
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  10. /*  $Id: blast_seqalign.cpp,v 1000.6 2004/06/01 18:05:52 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:  Christiam Camacho
  36. *
  37. * ===========================================================================
  38. */
  40. /// @file blast_seqalign.cpp
  41. /// Utility function to convert internal BLAST result structures into
  42. /// CSeq_align_set objects.
  44. #include <ncbi_pch.hpp>
  45. #include "blast_seqalign.hpp"
  47. #include <objects/seqloc/Seq_loc.hpp>
  48. #include <objects/seqloc/Seq_interval.hpp>
  49. #include <objects/seqalign/Seq_align.hpp>
  50. #include <objects/seqalign/Seq_align_set.hpp>
  51. #include <objects/seqalign/Dense_seg.hpp>
  52. #include <objects/seqalign/Dense_diag.hpp>
  53. #include <objects/seqalign/Std_seg.hpp>
  54. #include <objects/seqalign/Score.hpp>
  56. #include <objects/general/Object_id.hpp>
  57. #include <objmgr/util/sequence.hpp>
  58. #include <serial/iterator.hpp>
  60. /** @addtogroup AlgoBlast
  61.  *
  62.  * @{
  63.  */
  65. BEGIN_NCBI_SCOPE
  66. USING_SCOPE(objects);
  67. BEGIN_SCOPE(blast)
  69. #ifndef SMALLEST_EVALUE
  70. #define SMALLEST_EVALUE 1.0e-180
  71. #endif
  72. #ifndef GAP_VALUE
  73. #define GAP_VALUE -1
  74. #endif
  76. // Converts a frame into the appropriate strand
  77. static ENa_strand
  78. x_Frame2Strand(short frame)
  79. {
  80.     if (frame > 0)
  81.         return eNa_strand_plus;
  82.     else if (frame < 0)
  83.         return eNa_strand_minus;
  84.     else
  85.         return eNa_strand_unknown;
  86. }
  88. static int 
  89. x_GetCurrPos(int& pos, int pos2advance)
  90. {
  91.     int retval;
  93.     if (pos < 0)
  94.         retval = -(pos + pos2advance - 1);
  95.     else
  96.         retval = pos;
  97.     pos += pos2advance;
  98.     return retval;
  99. }
  101. static TSeqPos
  102. x_GetAlignmentStart(int& curr_pos, const GapEditScript* esp, 
  103.         ENa_strand strand, bool translate, int length, int original_length, 
  104.         short frame)
  105. {
  106.     TSeqPos retval;
  108.     if (strand == eNa_strand_minus) {
  110.         if (translate)
  111.             retval = original_length - 
  112.                 CODON_LENGTH*(x_GetCurrPos(curr_pos, esp->num) + esp->num)
  113.                 + frame + 1;
  114.         else
  115.             retval = length - x_GetCurrPos(curr_pos, esp->num) - esp->num;
  117.     } else {
  119.         if (translate)
  120.             retval = frame - 1 + CODON_LENGTH*x_GetCurrPos(curr_pos, esp->num);
  121.         else
  122.             retval = x_GetCurrPos(curr_pos, esp->num);
  124.     }
  126.     return retval;
  127. }
  129. /// C++ version of GXECollectDataForSeqalign
  130. /// Note that even though the edit_block is passed in, data for seqalign is
  131. /// collected from the esp argument for nsegs segments
  132. static void
  133. x_CollectSeqAlignData(const GapEditBlock* edit_block, 
  134.         const GapEditScript* esp_head, unsigned int nsegs,
  135.         vector<TSignedSeqPos>& starts, vector<TSeqPos>& lengths, 
  136.         vector<ENa_strand>& strands)
  137. {
  138.     ASSERT(edit_block != NULL);
  140.     GapEditScript* esp = const_cast<GapEditScript*>(esp_head);
  141.     ENa_strand m_strand, s_strand;      // strands of alignment
  142.     TSignedSeqPos m_start, s_start;     // running starts of alignment
  143.     int start1 = edit_block->start1;    // start of alignment on master sequence
  144.     int start2 = edit_block->start2;    // start of alignment on slave sequence
  146.     m_strand = x_Frame2Strand(edit_block->frame1);
  147.     s_strand = x_Frame2Strand(edit_block->frame2);
  149.     for (unsigned int i = 0; esp && i < nsegs; esp = esp->next, i++) {
  150.         switch (esp->op_type) {
  151.         case eGapAlignDecline:
  152.         case eGapAlignSub:
  153.             m_start = x_GetAlignmentStart(start1, esp, m_strand,
  154.                     edit_block->translate1 != 0, edit_block->length1,
  155.                     edit_block->original_length1, edit_block->frame1);
  157.             s_start = x_GetAlignmentStart(start2, esp, s_strand,
  158.                     edit_block->translate2 != 0, edit_block->length2,
  159.                     edit_block->original_length2, edit_block->frame2);
  161.             if (edit_block->reverse) {
  162.                 strands.push_back(s_strand);
  163.                 strands.push_back(m_strand);
  164.                 starts.push_back(s_start);
  165.                 starts.push_back(m_start);
  166.             } else {
  167.                 strands.push_back(m_strand);
  168.                 strands.push_back(s_strand);
  169.                 starts.push_back(m_start);
  170.                 starts.push_back(s_start);
  171.             }
  172.             break;
  174.         // Insertion on the master sequence (gap on slave)
  175.         case eGapAlignIns:
  176.             m_start = x_GetAlignmentStart(start1, esp, m_strand,
  177.                     edit_block->translate1 != 0, edit_block->length1,
  178.                     edit_block->original_length1, edit_block->frame1);
  180.             s_start = GAP_VALUE;
  182.             if (edit_block->reverse) {
  183.                 strands.push_back(i == 0 ? eNa_strand_unknown : s_strand);
  184.                 strands.push_back(m_strand);
  185.                 starts.push_back(s_start);
  186.                 starts.push_back(m_start);
  187.             } else {
  188.                 strands.push_back(m_strand);
  189.                 strands.push_back(i == 0 ? eNa_strand_unknown : s_strand);
  190.                 starts.push_back(m_start);
  191.                 starts.push_back(s_start);
  192.             }
  193.             break;
  195.         // Deletion on master sequence (gap; insertion on slave)
  196.         case eGapAlignDel:
  197.             m_start = GAP_VALUE;
  199.             s_start = x_GetAlignmentStart(start2, esp, s_strand,
  200.                     edit_block->translate2 != 0, edit_block->length2,
  201.                     edit_block->original_length2, edit_block->frame2);
  203.             if (edit_block->reverse) {
  204.                 strands.push_back(s_strand);
  205.                 strands.push_back(i == 0 ? eNa_strand_unknown : m_strand);
  206.                 starts.push_back(s_start);
  207.                 starts.push_back(m_start);
  208.             } else {
  209.                 strands.push_back(i == 0 ? eNa_strand_unknown : m_strand);
  210.                 strands.push_back(s_strand);
  211.                 starts.push_back(m_start);
  212.                 starts.push_back(s_start);
  213.             }
  214.             break;
  216.         default:
  217.             break;
  218.         }
  220.         lengths.push_back(esp->num);
  221.     }
  223.     // Make sure the vectors have the right size
  224.     if (lengths.size() != nsegs)
  225.         lengths.resize(nsegs);
  227.     if (starts.size() != nsegs*2)
  228.         starts.resize(nsegs*2);
  230.     if (strands.size() != nsegs*2)
  231.         strands.resize(nsegs*2);
  232. }
  234. static CRef<CDense_seg>
  235. x_CreateDenseg(const CSeq_id* master, const CSeq_id* slave,
  236.         vector<TSignedSeqPos>& starts, vector<TSeqPos>& lengths, 
  237.         vector<ENa_strand>& strands)
  238. {
  239.     ASSERT(master);
  240.     ASSERT(slave);
  242.     CRef<CDense_seg> dense_seg(new CDense_seg());
  244.     // Pairwise alignment is 2 dimensional
  245.     dense_seg->SetDim(2);
  247.     // Set the sequence ids
  248.     CDense_seg::TIds& ids = dense_seg->SetIds();
  249.     CRef<CSeq_id> tmp(new CSeq_id(master->AsFastaString()));
  250.     ids.push_back(tmp);
  251.     tmp.Reset(new CSeq_id(slave->AsFastaString()));
  252.     ids.push_back(tmp);
  253.     ids.resize(dense_seg->GetDim());
  255.     dense_seg->SetLens() = lengths;
  256.     dense_seg->SetStrands() = strands;
  257.     dense_seg->SetStarts() = starts;
  258.     dense_seg->SetNumseg(lengths.size());
  260.     return dense_seg;
  261. }
  263. static CSeq_align::C_Segs::TStd
  264. x_CreateStdSegs(const CSeq_id* master, const CSeq_id* slave, 
  265.         vector<TSignedSeqPos>& starts, vector<TSeqPos>& lengths, 
  266.         vector<ENa_strand>& strands, bool reverse, bool translate_master, bool
  267.         translate_slave)
  268. {
  269.     ASSERT(master);
  270.     ASSERT(slave);
  272.     CSeq_align::C_Segs::TStd retval;
  273.     int nsegs = lengths.size();         // number of segments in alignment
  274.     TSignedSeqPos m_start, m_stop;      // start and stop for master sequence
  275.     TSignedSeqPos s_start, s_stop;      // start and stop for slave sequence
  277.     CRef<CSeq_id> master_id(new CSeq_id(master->AsFastaString()));
  278.     CRef<CSeq_id> slave_id(new CSeq_id(slave->AsFastaString()));
  279.     
  280.     for (int i = 0; i < nsegs; i++) {
  281.         CRef<CStd_seg> std_seg(new CStd_seg());
  282.         CRef<CSeq_loc> master_loc(new CSeq_loc());
  283.         CRef<CSeq_loc> slave_loc(new CSeq_loc());
  285.         // Pairwise alignment is 2 dimensional
  286.         std_seg->SetDim(2);
  288.         // Set master seqloc
  289.         if ( (m_start = starts[2*i]) != GAP_VALUE) {
  290.             master_loc->SetInt().SetId(*master_id);
  291.             master_loc->SetInt().SetFrom(m_start);
  292.             if (translate_master || (reverse && translate_slave))
  293.                 m_stop = m_start + CODON_LENGTH*lengths[i] - 1;
  294.             else
  295.                 m_stop = m_start + lengths[i] - 1;
  296.             master_loc->SetInt().SetTo(m_stop);
  297.             master_loc->SetInt().SetStrand(strands[2*i]);
  298.         } else {
  299.             master_loc->SetEmpty(*master_id);
  300.         }
  302.         // Set slave seqloc
  303.         if ( (s_start = starts[2*i+1]) != GAP_VALUE) {
  304.             slave_loc->SetInt().SetId(*slave_id);
  305.             slave_loc->SetInt().SetFrom(s_start);
  306.             if (translate_slave || (reverse && translate_master))
  307.                 s_stop = s_start + CODON_LENGTH*lengths[i] - 1;
  308.             else
  309.                 s_stop = s_start + lengths[i] - 1;
  310.             slave_loc->SetInt().SetTo(s_stop);
  311.             slave_loc->SetInt().SetStrand(strands[2*i+1]);
  312.         } else {
  313.             slave_loc->SetEmpty(*slave_id);
  314.         }
  316.         if (reverse) {
  317.             std_seg->SetIds().push_back(slave_id);
  318.             std_seg->SetIds().push_back(master_id);
  319.             std_seg->SetLoc().push_back(slave_loc);
  320.             std_seg->SetLoc().push_back(master_loc);
  321.         } else {
  322.             std_seg->SetIds().push_back(master_id);
  323.             std_seg->SetIds().push_back(slave_id);
  324.             std_seg->SetLoc().push_back(master_loc);
  325.             std_seg->SetLoc().push_back(slave_loc);
  326.         }
  328.         retval.push_back(std_seg);
  329.     }
  331.     return retval;
  332. }
  334. /// Clone of GXECorrectUASequence (tools/gapxdrop.c)
  335. /// Assumes eGapAlignDecline regions are to the right of eGapAlign{Ins,Del}.
  336. /// This function swaps them (eGapAlignDecline ends to the right of the gap)
  337. static void 
  338. x_CorrectUASequence(GapEditBlock* edit_block)
  339. {
  340.     GapEditScript* curr = NULL,* curr_last = NULL;
  341.     GapEditScript* indel_prev = NULL; // pointer to node before the last
  342.             // insertion or deletion followed immediately by GAPALIGN_DECLINE
  343.     bool last_indel = false;    // last operation was an insertion or deletion
  345.     for (curr = edit_block->esp; curr; curr = curr->next) {
  347.         // if GAPALIGN_DECLINE immediately follows an insertion or deletion
  348.         if (curr->op_type == eGapAlignDecline && last_indel) {
  349.             /* This is invalid condition and regions should be
  350.                exchanged */
  352.             if (indel_prev != NULL)
  353.                 indel_prev->next = curr;
  354.             else
  355.                 edit_block->esp = curr; /* First element in a list */
  357.             // CC: If flaking gaps are allowed, curr_last could be NULL and the
  358.             // following statement would core dump...
  359.             curr_last->next = curr->next;
  360.             curr->next = curr_last;
  361.         }
  363.         last_indel = false;
  365.         if (curr->op_type == eGapAlignIns || curr->op_type == eGapAlignDel) {
  366.             last_indel = true;
  367.             indel_prev = curr_last;
  368.         }
  370.         curr_last = curr;
  371.     }
  372.     return;
  373. }
  375. /// C++ version of GXEMakeSeqAlign (tools/gapxdrop.c)
  376. /// Creates a Seq-align for a single HSP
  377. static CRef<CSeq_align>
  378. x_CreateSeqAlign(const CSeq_id* master, const CSeq_id* slave,
  379.     vector<TSignedSeqPos> starts, vector<TSeqPos> lengths,
  380.     vector<ENa_strand> strands, bool translate_master, bool translate_slave,
  381.     bool reverse)
  382. {
  383.     CRef<CSeq_align> sar(new CSeq_align());
  384.     sar->SetType(CSeq_align::eType_partial);
  385.     sar->SetDim(2);         // BLAST only creates pairwise alignments
  387.     if (translate_master || translate_slave) {
  388.         sar->SetSegs().SetStd() =
  389.             x_CreateStdSegs(master, slave, starts, lengths, strands,
  390.                     reverse, translate_master, translate_slave);
  391.     } else {
  392.         CRef<CDense_seg> dense_seg = 
  393.             x_CreateDenseg(master, slave, starts, lengths, strands);
  394.         sar->SetSegs().SetDenseg(*dense_seg);
  395.     }
  397.     return sar;
  398. }
  400. static CRef<CSeq_align>
  401. x_GapEditBlock2SeqAlign(GapEditBlock* edit_block, 
  402.         const CSeq_id* id1, const CSeq_id* id2)
  403. {
  404.     ASSERT(edit_block != NULL);
  406.     vector<TSignedSeqPos> starts;
  407.     vector<TSeqPos> lengths;
  408.     vector<ENa_strand> strands;
  409.     bool is_disc_align = false;
  410.     int nsegs = 0;      // number of segments in edit_block->esp
  412.     for (GapEditScript* t = edit_block->esp; t; t = t->next, nsegs++) {
  413.         if (t->op_type == eGapAlignDecline)
  414.             is_disc_align = true;
  415.     }
  417.     if (is_disc_align) {
  419.         /* By request of Steven Altschul - we need to have 
  420.            the unaligned part being to the left if it is adjacent to the
  421.            gap (insertion or deletion) - so this function will do
  422.            shuffeling */
  423.         x_CorrectUASequence(edit_block);
  425.         CRef<CSeq_align> seqalign(new CSeq_align());
  426.         seqalign->SetType(CSeq_align::eType_partial);
  427.         seqalign->SetDim(2);         // BLAST only creates pairwise alignments
  429.         bool skip_region;
  430.         GapEditScript* curr = NULL,* curr_head = edit_block->esp;
  432.         while (curr_head) {
  433.             skip_region = false;
  435.             for (nsegs = 0, curr = curr_head; curr; curr = curr->next, nsegs++){
  436.                 if (curr->op_type == eGapAlignDecline) {
  437.                     if (nsegs != 0) { // end of aligned region
  438.                         break;
  439.                     } else {
  440.                         while (curr && curr->op_type == eGapAlignDecline) {
  441.                             nsegs++;
  442.                             curr = curr->next;
  443.                         }
  444.                         skip_region = true;
  445.                         break;
  446.                     }
  447.                 }
  448.             }
  450.             // build seqalign for required regions only
  451.             if (!skip_region) {
  453.                 x_CollectSeqAlignData(edit_block, curr_head, nsegs, starts,
  454.                         lengths, strands);
  456.                 CRef<CSeq_align> sa_tmp = x_CreateSeqAlign(id1, id2, starts,
  457.                         lengths, strands, edit_block->translate1 !=0,
  458.                         edit_block->translate2 != 0, edit_block->reverse != 0);
  460.                 // Add this seqalign to the list
  461.                 if (sa_tmp)
  462.                     seqalign->SetSegs().SetDisc().Set().push_back(sa_tmp);
  463.             }
  464.             curr_head = curr;
  465.         }
  467.         return seqalign;
  469.     } else {
  471.         x_CollectSeqAlignData(edit_block, edit_block->esp, nsegs, 
  472.                 starts, lengths, strands);
  474.         return x_CreateSeqAlign(id1, id2, starts, lengths, strands,
  475.                 edit_block->translate1 != 0, edit_block->translate2 != 0,
  476.                 edit_block->reverse != 0);
  477.     }
  478. }
  480. /** Get the current position. */
  481. static Int4 get_current_pos(Int4* pos, Int4 length)
  482. {
  483.     Int4 val;
  484.     if(*pos < 0)
  485.         val = -(*pos + length -1);
  486.     else
  487.         val = *pos;
  488.     *pos += length;
  489.     return val;
  490. }
  492. /** This function is used for Out-Of-Frame traceback conversion
  493.  * Convert an OOF EditScript chain to a SeqAlign of type StdSeg.
  494.  * Used for a non-simple interval (i.e., one without subs. or 
  495.  * deletions).  
  496.  * The first SeqIdPtr in the chain of subject_id and query_id is 
  497.  * duplicated for the SeqAlign.
  498. */
  499. static CRef<CSeq_align>
  500. x_OOFEditBlock2SeqAlign(EProgram program, 
  501.     GapEditBlock* edit_block, 
  502.     const CSeq_id* query_id, const CSeq_id* subject_id)
  503. {
  504.     ASSERT(edit_block != NULL);
  506.     CRef<CSeq_align> seqalign(new CSeq_align());
  508.     Boolean reverse = FALSE;
  509.     GapEditScript* curr,* esp;
  510.     Int2 frame1, frame2;
  511.     Int4 start1, start2;
  512.     Int4 original_length1, original_length2;
  513.     CRef<CSeq_interval> seq_int1_last;
  514.     CRef<CSeq_interval> seq_int2_last;
  515.     CConstRef<CSeq_id> id1;
  516.     CConstRef<CSeq_id> id2;
  517.     CRef<CSeq_loc> slp1, slp2;
  518.     ENa_strand strand1, strand2;
  519.     bool first_shift;
  520.     Int4 from1, from2, to1, to2;
  521.     CRef<CSeq_id> tmp;
  523.     if (program == eBlastx) {
  524.        reverse = TRUE;
  525.        start1 = edit_block->start2;
  526.        start2 = edit_block->start1;
  527.        frame1 = edit_block->frame2;
  528.        frame2 = edit_block->frame1;
  529.        original_length1 = edit_block->original_length2;
  530.        original_length2 = edit_block->original_length1;
  531.        id1.Reset(subject_id);
  532.        id2.Reset(query_id);
  533.     } else { 
  534.        start1 = edit_block->start1;
  535.        start2 = edit_block->start2;
  536.        frame1 = edit_block->frame1;
  537.        frame2 = edit_block->frame2;
  538.        original_length1 = edit_block->original_length1;
  539.        original_length2 = edit_block->original_length2;
  540.        id1.Reset(query_id);
  541.        id2.Reset(subject_id);
  542.     }
  543.  
  544.     strand1 = x_Frame2Strand(frame1);
  545.     strand2 = x_Frame2Strand(frame2);
  546.     
  547.     esp = edit_block->esp;
  548.     
  549.     seqalign->SetDim(2); /**only two dimention alignment**/
  550.     
  551.     seqalign->SetType(CSeq_align::eType_partial); /**partial for gapped translating search. */
  553.     esp = edit_block->esp;
  555.     first_shift = false;
  557.     for (curr=esp; curr; curr=curr->next) {
  559.         slp1.Reset(new CSeq_loc());
  560.         slp2.Reset(new CSeq_loc());
  561.         
  562.         switch (curr->op_type) {
  563.         case eGapAlignDel: /* deletion of three nucleotides. */
  564.             
  565.             first_shift = false;
  567.             slp1->SetInt().SetFrom(get_current_pos(&start1, curr->num));
  568.             slp1->SetInt().SetTo(MIN(start1,original_length1) - 1);
  569.             tmp.Reset(new CSeq_id(id1->AsFastaString()));
  570.             slp1->SetInt().SetId(*tmp);
  571.             slp1->SetInt().SetStrand(strand1);
  572.             
  573.             /* Empty nucleotide piece */
  574.             tmp.Reset(new CSeq_id(id2->AsFastaString()));
  575.             slp2->SetEmpty(*tmp);
  577.             seq_int1_last.Reset(&slp1->SetInt());
  578.             /* Keep previous seq_int2_last, in case there is a frame shift
  579.                immediately after this gap */
  580.             
  581.             break;
  583.         case eGapAlignIns: /* insertion of three nucleotides. */
  584.             /* If gap is followed after frameshift - we have to
  585.                add this element for the alignment to be correct */
  586.             
  587.             if(first_shift) { /* Second frameshift in a row */
  588.                 /* Protein coordinates */
  589.                 slp1->SetInt().SetFrom(get_current_pos(&start1, 1));
  590.                 to1 = MIN(start1,original_length1) - 1;
  591.                 slp1->SetInt().SetTo(to1);
  592.                 tmp.Reset(new CSeq_id(id1->AsFastaString()));
  593.                 slp1->SetInt().SetId(*tmp);
  594.                 slp1->SetInt().SetStrand(strand1);
  595.                 
  596.                 /* Nucleotide scale shifted by op_type */
  597.                 from2 = get_current_pos(&start2, 3);
  598.                 to2 = MIN(start2,original_length2) - 1;
  599.                 slp2->SetInt().SetFrom(from2);
  600.                 slp2->SetInt().SetTo(to2);
  601.                 if (start2 > original_length2)
  602.                     slp1->SetInt().SetTo(to1 - 1);
  603.                 
  604.                 /* Transfer to DNA minus strand coordinates */
  605.                 if(strand2 == eNa_strand_minus) {
  606.                     slp2->SetInt().SetTo(original_length2 - from2 - 1);
  607.                     slp2->SetInt().SetFrom(original_length2 - to2 - 1);
  608.                 }
  609.                 
  610.                 tmp.Reset(new CSeq_id(id2->AsFastaString()));
  611.                 slp2->SetInt().SetId(*tmp);
  612.                 slp2->SetInt().SetStrand(strand2);
  614.                 CRef<CStd_seg> seg(new CStd_seg());
  615.                 seg->SetDim(2);
  617.                 CStd_seg::TIds& ids = seg->SetIds();
  619.                 if (reverse) {
  620.                     seg->SetLoc().push_back(slp2);
  621.                     seg->SetLoc().push_back(slp1);
  622.                     tmp.Reset(new CSeq_id(id2->AsFastaString()));
  623.                     ids.push_back(tmp);
  624.                     tmp.Reset(new CSeq_id(id1->AsFastaString()));
  625.                     ids.push_back(tmp);
  626.                 } else {
  627.                     seg->SetLoc().push_back(slp1);
  628.                     seg->SetLoc().push_back(slp2);
  629.                     tmp.Reset(new CSeq_id(id1->AsFastaString()));
  630.                     ids.push_back(tmp);
  631.                     tmp.Reset(new CSeq_id(id2->AsFastaString()));
  632.                     ids.push_back(tmp);
  633.                 }
  634.                 ids.resize(seg->GetDim());
  635.                 
  636.                 seqalign->SetSegs().SetStd().push_back(seg);
  637.             }
  639.             first_shift = false;
  641.             /* Protein piece is empty */
  642.             tmp.Reset(new CSeq_id(id1->AsFastaString()));
  643.             slp1->SetEmpty(*tmp);
  644.             
  645.             /* Nucleotide scale shifted by 3, protein gapped */
  646.             from2 = get_current_pos(&start2, curr->num*3);
  647.             to2 = MIN(start2,original_length2) - 1;
  648.             slp2->SetInt().SetFrom(from2);
  649.             slp2->SetInt().SetTo(to2);
  651.             /* Transfer to DNA minus strand coordinates */
  652.             if(strand2 == eNa_strand_minus) {
  653.                 slp2->SetInt().SetTo(original_length2 - from2 - 1);
  654.                 slp2->SetInt().SetFrom(original_length2 - to2 - 1);
  655.             }
  656.             tmp.Reset(new CSeq_id(id2->AsFastaString()));
  657.             slp2->SetInt().SetId(*tmp);
  658.             slp2->SetInt().SetStrand(strand2);
  659.             
  660.             seq_int1_last.Reset(NULL);
  661.             seq_int2_last.Reset(&slp2->SetInt()); /* Will be used to adjust "to" value */
  662.             
  663.             break;
  665.         case eGapAlignSub: /* Substitution. */
  667.             first_shift = false;
  669.             /* Protein coordinates */
  670.             from1 = get_current_pos(&start1, curr->num);
  671.             to1 = MIN(start1, original_length1) - 1;
  672.             /* Adjusting last segment and new start point in
  673.                nucleotide coordinates */
  674.             from2 = get_current_pos(&start2, curr->num*((Uint1)curr->op_type));
  675.             to2 = start2 - 1;
  676.             /* Chop off three bases and one residue at a time.
  677.                Why does this happen, seems like a bug?
  678.             */
  679.             while (to2 >= original_length2) {
  680.                 to2 -= 3;
  681.                 to1--;
  682.             }
  683.             /* Transfer to DNA minus strand coordinates */
  684.             if(strand2 == eNa_strand_minus) {
  685.                 int tmp_int;
  686.                 tmp_int = to2;
  687.                 to2 = original_length2 - from2 - 1;
  688.                 from2 = original_length2 - tmp_int - 1;
  689.             }
  691.             slp1->SetInt().SetFrom(from1);
  692.             slp1->SetInt().SetTo(to1);
  693.             tmp.Reset(new CSeq_id(id1->AsFastaString()));
  694.             slp1->SetInt().SetId(*tmp);
  695.             slp1->SetInt().SetStrand(strand1);
  696.             slp2->SetInt().SetFrom(from2);
  697.             slp2->SetInt().SetTo(to2);
  698.             tmp.Reset(new CSeq_id(id2->AsFastaString()));
  699.             slp2->SetInt().SetId(*tmp);
  700.             slp2->SetInt().SetStrand(strand2);
  701.            
  703.             seq_int1_last.Reset(&slp1->SetInt()); /* Will be used to adjust "to" value */
  704.             seq_int2_last.Reset(&slp2->SetInt()); /* Will be used to adjust "to" value */
  705.             
  706.             break;
  707.         case eGapAlignDel2: /* gap of two nucleotides. */
  708.         case eGapAlignDel1: /* Gap of one nucleotide. */
  709.         case eGapAlignIns1: /* Insertion of one nucleotide. */
  710.         case eGapAlignIns2: /* Insertion of two nucleotides. */
  712.             if(first_shift) { /* Second frameshift in a row */
  713.                 /* Protein coordinates */
  714.                 from1 = get_current_pos(&start1, 1);
  715.                 to1 = MIN(start1,original_length1) - 1;
  717.                 /* Nucleotide scale shifted by op_type */
  718.                 from2 = get_current_pos(&start2, (Uint1)curr->op_type);
  719.                 to2 = start2 - 1;
  720.                 if (to2 >= original_length2) {
  721.                     to2 = original_length2 -1;
  722.                     to1--;
  723.                 }
  724.                 /* Transfer to DNA minus strand coordinates */
  725.                 if(strand2 == eNa_strand_minus) {
  726.                     int tmp_int;
  727.                     tmp_int = to2;
  728.                     to2 = original_length2 - from2 - 1;
  729.                     from2 = original_length2 - tmp_int - 1;
  730.                 }
  732.                 slp1->SetInt().SetFrom(from1);
  733.                 slp1->SetInt().SetTo(to1);
  734.                 tmp.Reset(new CSeq_id(id1->AsFastaString()));
  735.                 slp1->SetInt().SetId(*tmp);
  736.                 slp1->SetInt().SetStrand(strand1);
  737.                 slp2->SetInt().SetFrom(from2);
  738.                 slp2->SetInt().SetTo(to2);
  739.                 tmp.Reset(new CSeq_id(id2->AsFastaString()));
  740.                 slp2->SetInt().SetId(*tmp);
  741.                 slp2->SetInt().SetStrand(strand2);
  743.                 seq_int1_last.Reset(&slp1->SetInt()); 
  744.                 seq_int2_last.Reset(&slp2->SetInt()); 
  746.                 break;
  747.             }
  748.             
  749.             first_shift = true;
  751.             /* If this substitution is following simple frameshift
  752.                we do not need to start new segment, but may continue
  753.                old one */
  755.             if(seq_int2_last) {
  756.                 get_current_pos(&start2, curr->num*((Uint1)curr->op_type-3));
  757.                 if(strand2 != eNa_strand_minus) {
  758.                     seq_int2_last->SetTo(start2 - 1);
  759.                 } else {
  760.                     /* Transfer to DNA minus strand coordinates */
  761.                     seq_int2_last->SetFrom(original_length2 - start2);
  762.                 }
  764.                 /* Adjustment for multiple shifts - theoretically possible,
  765.                    but very improbable */
  766.                 if(seq_int2_last->GetFrom() > seq_int2_last->GetTo()) {
  767.                     
  768.                     if(strand2 != eNa_strand_minus) {
  769.                         seq_int2_last->SetTo(seq_int2_last->GetTo() + 3);
  770.                     } else {
  771.                         seq_int2_last->SetFrom(seq_int2_last->GetFrom() - 3);
  772.                     }
  773.                     
  774.                     if (seq_int1_last.GetPointer() &&
  775. seq_int1_last->GetTo() != 0)
  776.                         seq_int1_last->SetTo(seq_int1_last->GetTo() + 1);
  777.                 }
  779.             } else if ((Uint1)curr->op_type > 3) {
  780.                 /* Protein piece is empty */
  781.                 tmp.Reset(new CSeq_id(id1->AsFastaString()));
  782.                 slp1->SetEmpty(*tmp);
  783.                 /* Simulating insertion of nucleotides */
  784.                 from2 = get_current_pos(&start2, 
  785.                                         curr->num*((Uint1)curr->op_type-3));
  786.                 to2 = MIN(start2,original_length2) - 1;
  788.                 /* Transfer to DNA minus strand coordinates */
  789.                 if(strand2 == eNa_strand_minus) {
  790.                     int tmp_int;
  791.                     tmp_int = to2;
  792.                     to2 = original_length2 - from2 - 1;
  793.                     from2 = original_length2 - tmp_int - 1;
  794.                 }
  795.                 slp2->SetInt().SetFrom(from2);
  796.                 slp2->SetInt().SetTo(to2);
  797.                 
  798.                 tmp.Reset(new CSeq_id(id2->AsFastaString()));
  799.                 slp2->SetInt().SetId(*tmp);
  801.                 seq_int1_last.Reset(NULL);
  802.                 seq_int2_last.Reset(&slp2->SetInt()); /* Will be used to adjust "to" value */
  803.                 break;
  804.             } else {
  805.                 continue;       /* Main loop */
  806.             }
  807.             continue;       /* Main loop */
  808.             /* break; */
  809.         default:
  810.             continue;       /* Main loop */
  811.             /* break; */
  812.         } 
  814.         CRef<CStd_seg> seg(new CStd_seg());
  815.         seg->SetDim(2);
  816.         
  817.         CStd_seg::TIds& ids = seg->SetIds();
  819.         if (reverse) {
  820.             seg->SetLoc().push_back(slp2);
  821.             seg->SetLoc().push_back(slp1);
  822.             tmp.Reset(new CSeq_id(id2->AsFastaString()));
  823.             ids.push_back(tmp);
  824.             tmp.Reset(new CSeq_id(id1->AsFastaString()));
  825.             ids.push_back(tmp);
  826.         } else {
  827.             seg->SetLoc().push_back(slp1);
  828.             seg->SetLoc().push_back(slp2);
  829.             tmp.Reset(new CSeq_id(id1->AsFastaString()));
  830.             ids.push_back(tmp);
  831.             tmp.Reset(new CSeq_id(id2->AsFastaString()));
  832.             ids.push_back(tmp);
  833.         }
  834.         ids.resize(seg->GetDim());
  835.         
  836.         seqalign->SetSegs().SetStd().push_back(seg);
  837.     }
  838.     
  839.     return seqalign;
  840. }
  842. /// Creates and initializes CScore with i (if it's non-zero) or d
  843. static CRef<CScore> 
  844. x_MakeScore(const string& ident_string, double d = 0.0, int i = 0)
  845. {
  846.     CRef<CScore> retval(new CScore());
  848.     CRef<CObject_id> id(new CObject_id());
  849.     id->SetStr(ident_string);
  850.     retval->SetId(*id);
  852.     CRef<CScore::C_Value> val(new CScore::C_Value());
  853.     if (i)
  854.         val->SetInt(i);
  855.     else
  856.         val->SetReal(d);
  857.     retval->SetValue(*val);
  859.     return retval;
  860. }
  862. /// C++ version of GetScoreSetFromBlastHsp (tools/blastutl.c)
  863. static void
  864. x_BuildScoreList(const BlastHSP* hsp, const BlastScoreBlk* sbp, const
  865.         BlastScoringOptions* score_options, CSeq_align::TScore& scores,
  866.         EProgram program)
  867. {
  868.     string score_type;
  869.     Blast_KarlinBlk* kbp = NULL;
  871.     if (!hsp)
  872.         return;
  874.     score_type = "score";
  875.     if (hsp->score)
  876.         scores.push_back(x_MakeScore(score_type, 0.0, hsp->score));
  878.     score_type = "sum_n";
  879.     if (hsp->num > 1)
  880.         scores.push_back(x_MakeScore(score_type, 0.0, hsp->num));
  882.     // Set the E-Value
  883.     double evalue = (hsp->evalue < SMALLEST_EVALUE) ? 0.0 : hsp->evalue;
  884.     score_type = (hsp->num == 1) ? "e_value" : "sum_e";
  885.     if (evalue >= 0.0)
  886.         scores.push_back(x_MakeScore(score_type, evalue));
  888.     // Calculate the bit score from the raw score
  889.     score_type = "bit_score";
  890.     if (program == eBlastn || !score_options->gapped_calculation)
  891.         kbp = sbp->kbp[hsp->context];
  892.     else
  893.         kbp = sbp->kbp_gap[hsp->context];
  895.     double bit_score = ((hsp->score*kbp->Lambda) - kbp->logK)/NCBIMATH_LN2;
  896.     if (bit_score >= 0.0)
  897.         scores.push_back(x_MakeScore(score_type, bit_score));
  899.     // Set the identity score
  900.     score_type = "num_ident";
  901.     if (hsp->num_ident > 0)
  902.         scores.push_back(x_MakeScore(score_type, 0.0, hsp->num_ident));
  904.     if (hsp->splice_junction > 0) {
  905.         // Splice junction(s) was (were) found between linked HSPs
  906.         score_type = "splice_junction";
  907.         scores.push_back(x_MakeScore(score_type, 0.0, hsp->splice_junction));
  908.     }
  910.     return;
  911. }
  914. static void
  915. x_AddScoresToSeqAlign(CRef<CSeq_align>& seqalign, const BlastHSP* hsp, 
  916.         const BlastScoreBlk* sbp, EProgram program,
  917.         const BlastScoringOptions* score_options)
  918. {
  919.     // Add the scores for this HSP
  920.     CSeq_align::TScore& score_list = seqalign->SetScore();
  921.     x_BuildScoreList(hsp, sbp, score_options, score_list, program);
  922. }
  925. CRef<CDense_diag>
  926. x_UngappedHSPToDenseDiag(BlastHSP* hsp, const CSeq_id *query_id, 
  927.     const CSeq_id *subject_id, const BlastScoreBlk* sbp,
  928.     const BlastScoringOptions* score_options, EProgram program, 
  929.     Int4 query_length, Int4 subject_length)
  930. {
  931.     CRef<CDense_diag> retval(new CDense_diag());
  932.     
  933.     // Pairwise alignment is 2 dimensional
  934.     retval->SetDim(2);
  936.     // Set the sequence ids
  937.     CDense_diag::TIds& ids = retval->SetIds();
  938.     CRef<CSeq_id> tmp(new CSeq_id(query_id->AsFastaString()));
  939.     ids.push_back(tmp);
  940.     tmp.Reset(new CSeq_id(subject_id->AsFastaString()));
  941.     ids.push_back(tmp);
  942.     ids.resize(retval->GetDim());
  944.     retval->SetLen(hsp->query.length);
  946.     CDense_diag::TStrands& strands = retval->SetStrands();
  947.     strands.push_back(x_Frame2Strand(hsp->query.frame));
  948.     strands.push_back(x_Frame2Strand(hsp->subject.frame));
  949.     CDense_diag::TStarts& starts = retval->SetStarts();
  950.     if (hsp->query.frame >= 0) {
  951.        starts.push_back(hsp->query.offset);
  952.     } else {
  953.        starts.push_back(query_length - hsp->query.offset - hsp->query.length);
  954.     }
  955.     if (hsp->subject.frame >= 0) {
  956.        starts.push_back(hsp->subject.offset);
  957.     } else {
  958.        starts.push_back(subject_length - hsp->subject.offset - 
  959.                         hsp->subject.length);
  960.     }
  962.     CSeq_align::TScore& score_list = retval->SetScores();
  963.     x_BuildScoreList(hsp, sbp, score_options, score_list, program);
  965.     return retval;
  966. }
  968. CRef<CStd_seg>
  969. x_UngappedHSPToStdSeg(BlastHSP* hsp, const CSeq_id *query_id, 
  970.     const CSeq_id *subject_id, const BlastScoreBlk* sbp,
  971.     const BlastScoringOptions* score_options, EProgram program, 
  972.     Int4 query_length, Int4 subject_length)
  973. {
  974.     CRef<CStd_seg> retval(new CStd_seg());
  976.     // Pairwise alignment is 2 dimensional
  977.     retval->SetDim(2);
  979.     CRef<CSeq_loc> query_loc(new CSeq_loc());
  980.     CRef<CSeq_loc> subject_loc(new CSeq_loc());
  982.     // Set the sequence ids
  983.     CStd_seg::TIds& ids = retval->SetIds();
  984.     CRef<CSeq_id> tmp(new CSeq_id(query_id->AsFastaString()));
  985.     query_loc->SetInt().SetId(*tmp);
  986.     ids.push_back(tmp);
  987.     tmp.Reset(new CSeq_id(subject_id->AsFastaString()));
  988.     subject_loc->SetInt().SetId(*tmp);
  989.     ids.push_back(tmp);
  990.     ids.resize(retval->GetDim());
  992.     query_loc->SetInt().SetStrand(x_Frame2Strand(hsp->query.frame));
  993.     subject_loc->SetInt().SetStrand(x_Frame2Strand(hsp->subject.frame));
  995.     if (hsp->query.frame == 0) {
  996.        query_loc->SetInt().SetFrom(hsp->query.offset);
  997.        query_loc->SetInt().SetTo(hsp->query.offset + hsp->query.length - 1);
  998.     } else if (hsp->query.frame > 0) {
  999.        query_loc->SetInt().SetFrom(CODON_LENGTH*(hsp->query.offset) + 
  1000.                                    hsp->query.frame - 1);
  1001.        query_loc->SetInt().SetTo(CODON_LENGTH*(hsp->query.offset+
  1002.                                                hsp->query.length)
  1003.                                  + hsp->query.frame - 2);
  1004.     } else {
  1005.        query_loc->SetInt().SetFrom(query_length -
  1006.            CODON_LENGTH*(hsp->query.offset+hsp->query.length) +
  1007.            hsp->query.frame + 1);
  1008.        query_loc->SetInt().SetTo(query_length - CODON_LENGTH*hsp->query.offset
  1009.                                  + hsp->query.frame);
  1010.     }
  1012.     if (hsp->subject.frame == 0) {
  1013.        subject_loc->SetInt().SetFrom(hsp->subject.offset);
  1014.        subject_loc->SetInt().SetTo(hsp->subject.offset + 
  1015.                                    hsp->subject.length - 1);
  1016.     } else if (hsp->subject.frame > 0) {
  1017.        subject_loc->SetInt().SetFrom(CODON_LENGTH*(hsp->subject.offset) + 
  1018.                                    hsp->subject.frame - 1);
  1019.        subject_loc->SetInt().SetTo(CODON_LENGTH*(hsp->subject.offset+
  1020.                                                hsp->subject.length) +
  1021.                                    hsp->subject.frame - 2);
  1023.     } else {
  1024.        subject_loc->SetInt().SetFrom(subject_length -
  1025.            CODON_LENGTH*(hsp->subject.offset+hsp->subject.length) +
  1026.            hsp->subject.frame + 1);
  1027.        subject_loc->SetInt().SetTo(subject_length - 
  1028.            CODON_LENGTH*hsp->subject.offset + hsp->subject.frame);
  1029.     }
  1031.     retval->SetLoc().push_back(query_loc);
  1032.     retval->SetLoc().push_back(subject_loc);
  1034.     CSeq_align::TScore& score_list = retval->SetScores();
  1035.     x_BuildScoreList(hsp, sbp, score_options, score_list, program);
  1037.     return retval;
  1038. }
  1040. static CRef<CSeq_align>
  1041. BLASTUngappedHspListToSeqAlign(EProgram program, 
  1042.     BlastHSPList* hsp_list, const CSeq_id *query_id, 
  1043.     const CSeq_id *subject_id, Int4 query_length, Int4 subject_length,
  1044.     const BlastScoringOptions* score_options, const BlastScoreBlk* sbp)
  1045. {
  1046.     CRef<CSeq_align> retval(new CSeq_align()); 
  1047.     BlastHSP** hsp_array;
  1048.     int index;
  1050.     retval->SetType(CSeq_align::eType_diags);
  1052.     hsp_array = hsp_list->hsp_array;
  1054.     /* All HSPs are put in one seqalign, containing a list of 
  1055.      * DenseDiag for same molecule search, or StdSeg for translated searches.
  1056.     */
  1057.     if (program == eBlastn || 
  1058.         program == eBlastp ||
  1059.         program == eRPSBlast) {
  1060.         for (index=0; index<hsp_list->hspcnt; index++) { 
  1061.             BlastHSP* hsp = hsp_array[index];
  1062.             retval->SetSegs().SetDendiag().push_back(
  1063.                 x_UngappedHSPToDenseDiag(hsp, query_id, subject_id, sbp, 
  1064.                     score_options, program, query_length, subject_length));
  1065.         }
  1066.     } else { // Translated search
  1067.         for (index=0; index<hsp_list->hspcnt; index++) { 
  1068.             BlastHSP* hsp = hsp_array[index];
  1069.             retval->SetSegs().SetStd().push_back(
  1070.                 x_UngappedHSPToStdSeg(hsp, query_id, subject_id, sbp, 
  1071.                     score_options, program, query_length, subject_length));
  1072.         }
  1073.     }
  1075.     return retval;
  1076. }
  1078. // This is called for each query and each subject in the BLAST search
  1079. // We always return CSeq_aligns of type disc to allow multiple HSPs
  1080. // corresponding to the same query-subject pair to be grouped in one CSeq_align
  1081. static CRef<CSeq_align>
  1082. BLASTHspListToSeqAlign(EProgram program, 
  1083.     BlastHSPList* hsp_list, const CSeq_id *query_id, 
  1084.     const CSeq_id *subject_id,
  1085.     const BlastScoringOptions* score_options, const BlastScoreBlk* sbp)
  1086. {
  1087.     CRef<CSeq_align> retval(new CSeq_align()); 
  1088.     retval->SetType(CSeq_align::eType_disc);
  1089.     retval->SetDim(2);         // BLAST only creates pairwise alignments
  1091.     // Process the list of HSPs corresponding to one subject sequence and
  1092.     // create one seq-align for each list of HSPs (use disc seqalign when
  1093.     // multiple HSPs are found).
  1094.     BlastHSP** hsp_array = hsp_list->hsp_array;
  1096.     for (int index = 0; index < hsp_list->hspcnt; index++) { 
  1097.         BlastHSP* hsp = hsp_array[index];
  1098.         CRef<CSeq_align> seqalign;
  1100.         if (score_options->is_ooframe) {
  1101.             seqalign = 
  1102.                 x_OOFEditBlock2SeqAlign(program, hsp->gap_info, 
  1103.                     query_id, subject_id);
  1104.         } else {
  1105.             seqalign = 
  1106.                 x_GapEditBlock2SeqAlign(hsp->gap_info, 
  1107.                     query_id, subject_id);
  1108.         }
  1109.         
  1110.         x_AddScoresToSeqAlign(seqalign, hsp, sbp, program, 
  1111.                               score_options);
  1112.         retval->SetSegs().SetDisc().Set().push_back(seqalign);
  1113.     }
  1114.     
  1115.     return retval;
  1116. }
  1118. static void
  1119. x_GetSequenceLengthAndId(const SSeqLoc* ss,          // [in]
  1120.                          const BlastSeqSrc* seq_src, // [in]
  1121.                          int oid,                    // [in] 
  1122.                          CSeq_id** seqid,            // [out]
  1123.                          TSeqPos* length)            // [out]
  1124. {
  1125.     ASSERT(ss || seq_src);
  1126.     ASSERT(seqid);
  1127.     ASSERT(length);
  1129.     if ( !seq_src ) {
  1130.         *seqid = const_cast<CSeq_id*>(&sequence::GetId(*ss->seqloc,
  1131.                                                        ss->scope));
  1132.         *length = sequence::GetLength(*ss->seqloc, ss->scope);
  1133.     } else {
  1134.         ListNode* seqid_wrap;
  1135.         seqid_wrap = BLASTSeqSrcGetSeqId(seq_src, (void*) &oid);
  1136.         ASSERT(seqid_wrap);
  1137.         if (seqid_wrap->choice == BLAST_SEQSRC_CPP_SEQID) {
  1138.             *seqid = (CSeq_id*)seqid_wrap->ptr;
  1139.             ListNodeFree(seqid_wrap);
  1140.         } else if (seqid_wrap->choice == BLAST_SEQSRC_CPP_SEQID_REF) {
  1141.             *seqid = ((CRef<CSeq_id>*)seqid_wrap->ptr)->GetPointer();
  1142.         } else {
  1143.             /** FIXME!!! This is wrong, because the id created here will 
  1144.                 not be registered! However if sequence source returns a 
  1145.                 C object, we cannot handle it here. */
  1146.             char* id = BLASTSeqSrcGetSeqIdStr(seq_src, (void*) &oid);
  1147.             string id_str(id);
  1148.             *seqid = new CSeq_id(id_str);
  1149.             sfree(id);
  1150.         }
  1151.         *length = BLASTSeqSrcGetSeqLen(seq_src, (void*) &oid);
  1152.     }
  1153.     return;
  1154. }
  1156. /// Constructs an empty seq-align-set containing an empty discontinuous
  1157. /// seq-align.
  1158. static CSeq_align_set*
  1159. x_CreateEmptySeq_align_set(CSeq_align_set* sas)
  1160. {
  1161.     CSeq_align_set* retval = NULL;
  1163.     if (!sas) {
  1164.         retval = new CSeq_align_set;
  1165.     } else {
  1166.         retval = sas;
  1167.     }
  1169.     CRef<CSeq_align> empty_seqalign(new CSeq_align);
  1170.     empty_seqalign->SetType(CSeq_align::eType_disc);
  1171.     empty_seqalign->SetSegs().SetDisc(*new CSeq_align_set);
  1173.     retval->Set().push_back(empty_seqalign);
  1174.     return retval;
  1175. }
  1177. /// Always remap the query, the subject is remapped if it's given (assumes
  1178. /// alignment created by BLAST 2 Sequences API).
  1179. /// Since the query strands were already taken into account when CSeq_align 
  1180. /// was created, only start position shifts in the CSeq_loc's are relevant in 
  1181. /// this function. However full remapping is necessary for the subject sequence
  1182. /// if it is on a negative strand.
  1184. static void
  1185. x_RemapAlignmentCoordinates(CRef<CSeq_align> sar, 
  1186.                             const SSeqLoc* query,
  1187.                             const SSeqLoc* subject = NULL)
  1188. {
  1189.     _ASSERT(sar);
  1190.     ASSERT(query);
  1191.     const int query_dimension = 0;
  1192.     const int subject_dimension = 1;
  1194.     // If subject is on a minus strand, we'll need to flip subject strands 
  1195.     // and remap subject coordinates on all segments.
  1196.     // Otherwise we only need to shift query and/or subject coordinates, 
  1197.     // if the respective location starts not from 0.
  1198.     bool remap_subject = 
  1199.         (subject && subject->seqloc->IsInt() &&
  1200.          subject->seqloc->GetInt().GetStrand() == eNa_strand_minus);
  1202.     TSeqPos q_shift = 0, s_shift = 0;
  1204.     if (query->seqloc->IsInt()) {
  1205.         q_shift = query->seqloc->GetInt().GetFrom();
  1206.     }
  1207.     if (subject && subject->seqloc->IsInt()) {
  1208.         s_shift = subject->seqloc->GetInt().GetFrom();
  1209.     }
  1211.     if (remap_subject || q_shift > 0 || s_shift > 0) {
  1212.         for (CTypeIterator<CDense_seg> itr(Begin(*sar)); itr; ++itr) {
  1213.             const vector<ENa_strand> strands = itr->GetStrands();
  1214.             // Create temporary CSeq_locs with strands either matching 
  1215.             // (for query and for subject if it is not on a minus strand),
  1216.             // or opposite to those in the segment, to force RemapToLoc to 
  1217.             // behave in the correct way.
  1218.             if (q_shift > 0) {
  1219.                 CSeq_loc q_seqloc;
  1220.                 ENa_strand q_strand = strands[0];
  1221.                 q_seqloc.SetInt().SetFrom(q_shift);
  1222.                 q_seqloc.SetInt().SetTo(query->seqloc->GetInt().GetTo());
  1223.                 q_seqloc.SetInt().SetStrand(q_strand);
  1224.                 q_seqloc.SetInt().SetId().Assign(sequence::GetId(*query->seqloc, query->scope));
  1225.                 itr->RemapToLoc(query_dimension, q_seqloc);
  1226.             }
  1227.             if (remap_subject || s_shift > 0) {
  1228.                 CSeq_loc s_seqloc;
  1229.                 ENa_strand s_strand;
  1230.                 if (remap_subject) {
  1231.                     s_strand = ((strands[1] == eNa_strand_plus) ? 
  1232.                                 eNa_strand_minus : eNa_strand_plus);
  1233.                 } else {
  1234.                     s_strand = strands[1];
  1235.                 }
  1236.                 s_seqloc.SetInt().SetFrom(s_shift);
  1237.                 s_seqloc.SetInt().SetTo(subject->seqloc->GetInt().GetTo());
  1238.                 s_seqloc.SetInt().SetStrand(s_strand);
  1239.                 s_seqloc.SetInt().SetId().Assign(sequence::GetId(*subject->seqloc, subject->scope));
  1240.                 itr->RemapToLoc(subject_dimension, s_seqloc);
  1241.             }
  1242.         }
  1243.     }
  1244. }
  1246. CSeq_align_set*
  1247. BLAST_HitList2CSeqAlign(const BlastHitList* hit_list, 
  1248.     EProgram prog, SSeqLoc &query,
  1249.     const BlastSeqSrc* seq_src,
  1250.     const BlastScoringOptions* score_options, 
  1251.     const BlastScoreBlk* sbp)
  1252. {
  1253.     CSeq_align_set* seq_aligns = new CSeq_align_set();
  1254.     bool is_gapped = score_options->gapped_calculation ? true : false;
  1256.     ASSERT(seq_src);
  1258.     if (!hit_list) {
  1259.         return x_CreateEmptySeq_align_set(seq_aligns);
  1260.     }
  1262.     TSeqPos query_length = 0;
  1263.     CSeq_id* qid = NULL;
  1264.     CConstRef<CSeq_id> query_id;
  1265.     x_GetSequenceLengthAndId(&query, NULL, 0, &qid, &query_length);
  1266.     query_id.Reset(qid);
  1268.     TSeqPos subj_length = 0;
  1269.     CSeq_id* sid = NULL;
  1270.     CConstRef<CSeq_id> subject_id;
  1272.     for (int index = 0; index < hit_list->hsplist_count; index++) {
  1273.         BlastHSPList* hsp_list = hit_list->hsplist_array[index];
  1274.         if (!hsp_list)
  1275.             continue;
  1277.         x_GetSequenceLengthAndId(NULL, seq_src, hsp_list->oid, 
  1278.                                  &sid, &subj_length);
  1279.         subject_id.Reset(sid);
  1280.         
  1281.         // Create a CSeq_align for each matching sequence
  1282.         CRef<CSeq_align> hit_align;
  1283.         if (is_gapped) {
  1284.             hit_align = 
  1285.                 BLASTHspListToSeqAlign(prog, hsp_list, query_id, 
  1286.                                        subject_id, score_options, sbp);
  1287.         } else {
  1288.             hit_align = 
  1289.                 BLASTUngappedHspListToSeqAlign(prog, hsp_list, query_id, 
  1290.                     subject_id, query_length, subj_length, score_options, sbp);
  1291.         }
  1292.         ListNode* subject_loc_wrap = 
  1293.             BLASTSeqSrcGetSeqLoc(seq_src, (void*)&hsp_list->oid);
  1294.         SSeqLoc* subject_loc = NULL;
  1295.         if (subject_loc_wrap && 
  1296.             subject_loc_wrap->choice == BLAST_SEQSRC_CPP_SEQLOC)
  1297.             subject_loc = (SSeqLoc*) subject_loc_wrap->ptr;
  1298.         x_RemapAlignmentCoordinates(hit_align, &query, subject_loc);
  1299.         seq_aligns->Set().push_back(hit_align);
  1300.     }
  1301.     return seq_aligns;
  1302. }
  1304. TSeqAlignVector
  1305. BLAST_Results2CSeqAlign(const BlastHSPResults* results, 
  1306.         EProgram prog,
  1307.         TSeqLocVector &query,
  1308.         const BlastSeqSrc* seq_src,
  1309.         const BlastScoringOptions* score_options, 
  1310.         const BlastScoreBlk* sbp)
  1311. {
  1312.     ASSERT(results->num_queries == (int)query.size());
  1313.     ASSERT(seq_src);
  1315.     TSeqAlignVector retval;
  1316.     CConstRef<CSeq_id> query_id;
  1318.     // Process each query's hit list
  1319.     for (int index = 0; index < results->num_queries; index++) {
  1320.        BlastHitList* hit_list = results->hitlist_array[index];
  1322.        CRef<CSeq_align_set> seq_aligns(BLAST_HitList2CSeqAlign(hit_list, prog,
  1323.                                            query[index], seq_src,
  1324.                                            score_options, sbp));
  1326.        retval.push_back(seq_aligns);
  1327.        _TRACE("Query " << index << ": " << seq_aligns->Get().size() 
  1328.               << " seqaligns");
  1330.     }
  1331.     
  1332.     return retval;
  1333. }
  1335. TSeqAlignVector
  1336. BLAST_OneSubjectResults2CSeqAlign(const BlastHSPResults* results, 
  1337.         EProgram prog,
  1338.         TSeqLocVector &query,
  1339.         const BlastSeqSrc* seq_src,
  1340.         Int4 subject_index,
  1341.         const BlastScoringOptions* score_options, 
  1342.         const BlastScoreBlk* sbp)
  1343. {
  1344.     ASSERT(results->num_queries == (int)query.size());
  1345.     ASSERT(seq_src);
  1347.     TSeqAlignVector retval;
  1348.     CConstRef<CSeq_id> subject_id;
  1349.     CConstRef<CSeq_id> query_id;
  1350.     CSeq_id* sid = NULL;
  1351.     CSeq_id* qid = NULL;
  1352.     bool is_gapped = score_options->gapped_calculation ? true : false;
  1353.     TSeqPos subj_length = 0;
  1354.     TSeqPos query_length = 0;
  1356.     // Subject is the same for all queries, so retrieve its id right away
  1357.     x_GetSequenceLengthAndId(NULL, seq_src, subject_index, 
  1358.                              &sid, &subj_length);
  1359.     subject_id.Reset(sid);
  1360.         
  1362.     // Process each query's hit list
  1363.     for (int index = 0; index < results->num_queries; index++) {
  1364.         x_GetSequenceLengthAndId(&query[index], NULL, 0, &qid, &query_length);
  1365.         query_id.Reset(qid);
  1366.         CRef<CSeq_align_set> seq_aligns;
  1367.         BlastHitList* hit_list = results->hitlist_array[index];
  1368.         BlastHSPList* hsp_list = NULL;
  1369.         // Find the HSP list corresponding to this subject, if it exists
  1370.         if (hit_list) {
  1371.             int result_index;
  1372.             for (result_index = 0; result_index < hit_list->hsplist_count;
  1373.                  ++result_index) {
  1374.                 hsp_list = hit_list->hsplist_array[result_index];
  1375.                 if (hsp_list->oid == subject_index)
  1376.                     break;
  1377.             }
  1378.         }
  1380.         if (hsp_list) {
  1381.             CRef<CSeq_align> hit_align;
  1382.             if (is_gapped) {
  1383.                 hit_align = 
  1384.                     BLASTHspListToSeqAlign(prog, hsp_list, query_id, 
  1385.                                            subject_id, score_options, sbp);
  1386.             } else {
  1387.                 hit_align = 
  1388.                     BLASTUngappedHspListToSeqAlign(prog, hsp_list, query_id, 
  1389.                         subject_id, query_length, subj_length, 
  1390.                         score_options, sbp);
  1391.             }
  1392.             ListNode* subject_loc_wrap = 
  1393.                 BLASTSeqSrcGetSeqLoc(seq_src, (void*)&hsp_list->oid);
  1394.             SSeqLoc* subject_loc = NULL;
  1395.             if (subject_loc_wrap && 
  1396.                 subject_loc_wrap->choice == BLAST_SEQSRC_CPP_SEQLOC)
  1397.                 subject_loc = (SSeqLoc*) subject_loc_wrap->ptr;
  1398.             x_RemapAlignmentCoordinates(hit_align, &query[index], 
  1399.                                         subject_loc);
  1400.             ListNodeFree(subject_loc_wrap);
  1401.             seq_aligns.Reset(new CSeq_align_set());
  1402.             seq_aligns->Set().push_back(hit_align);
  1403.         } else {
  1404.             seq_aligns.Reset(x_CreateEmptySeq_align_set(NULL));
  1405.         }
  1406.         retval.push_back(seq_aligns);
  1407.     }
  1408.     
  1409.     return retval;
  1410. }
  1412. END_SCOPE(blast)
  1413. END_NCBI_SCOPE
  1415. /* @} */
  1417. /*
  1418. * ===========================================================================
  1419. *
  1420. * $Log: blast_seqalign.cpp,v $
  1421. * Revision 1000.6  2004/06/01 18:05:52  gouriano
  1422. * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.42
  1423. *
  1424. * Revision 1.42  2004/05/21 21:41:02  gorelenk
  1425. * Added PCH ncbi_pch.hpp
  1426. *
  1427. * Revision 1.41  2004/05/19 15:26:04  dondosha
  1428. * Edit script operations changed from macros to an enum
  1429. *
  1430. * Revision 1.40  2004/05/05 14:42:25  dondosha
  1431. * Correction in x_RemapAlignmentCoordinates for whole vs. interval Seq-locs
  1432. *
  1433. * Revision 1.39  2004/04/28 19:40:45  dondosha
  1434. * Fixed x_RemapAlignmentCoordinates function to work correctly with all strand combinations
  1435. *
  1436. * Revision 1.38  2004/04/19 12:59:12  madden
  1437. * Changed BLAST_KarlinBlk to Blast_KarlinBlk to avoid conflict with blastkar.h structure
  1438. *
  1439. * Revision 1.37  2004/04/16 14:28:19  papadopo
  1440. * add use of eRPSBlast program
  1441. *
  1442. * Revision 1.36  2004/04/06 20:47:14  dondosha
  1443. * Check if BLASTSeqSrcGetSeqId returns a pointer to CRef instead of a simple pointer to CSeq_id
  1444. *
  1445. * Revision 1.35  2004/03/24 19:14:14  dondosha
  1446. * BlastHSP structure does not have ordering_method field any more, but it does contain a splice_junction field
  1447. *
  1448. * Revision 1.34  2004/03/23 18:22:06  dondosha
  1449. * Minor memory leak fix
  1450. *
  1451. * Revision 1.33  2004/03/23 14:10:50  camacho
  1452. * Minor doxygen fix
  1453. *
  1454. * Revision 1.32  2004/03/19 19:22:55  camacho
  1455. * Move to doxygen group AlgoBlast, add missing CVS logs at EOF
  1456. *
  1457. * Revision 1.31  2004/03/16 22:03:38  camacho
  1458. * Remove dead code
  1459. *
  1460. * Revision 1.30  2004/03/15 19:58:55  dondosha
  1461. * Added BLAST_OneSubjectResults2CSeqalign function to retrieve single subject results from BlastHSPResults
  1462. *
  1463. * Revision 1.29  2003/12/19 20:16:10  dondosha
  1464. * Get length in x_GetSequenceLengthAndId regardless of whether search is gapped; do not call RemapToLoc for whole Seq-locs
  1465. *
  1466. * Revision 1.28  2003/12/03 16:43:47  dondosha
  1467. * Renamed BlastMask to BlastMaskLoc, BlastResults to BlastHSPResults
  1468. *
  1469. * Revision 1.27  2003/12/01 20:02:39  coulouri
  1470. * fix msvc warning
  1471. *
  1472. * Revision 1.26  2003/11/24 20:59:12  ucko
  1473. * Change one ASSERT to _ASSERT (probably more appropriate in general)
  1474. * due to MSVC brokenness.
  1475. *
  1476. * Revision 1.25  2003/11/24 17:14:58  camacho
  1477. * Remap alignment coordinates to original Seq-locs
  1478. *
  1479. * Revision 1.24  2003/11/04 18:37:36  dicuccio
  1480. * Fix for brain-dead MSVC (operator && is ambiguous)
  1481. *
  1482. * Revision 1.23  2003/11/04 17:13:31  dondosha
  1483. * Implemented conversion of results to seqalign for out-of-frame search
  1484. *
  1485. * Revision 1.22  2003/10/31 22:08:39  dondosha
  1486. * Implemented conversion of BLAST results to seqalign for ungapped search
  1487. *
  1488. * Revision 1.21  2003/10/31 00:05:15  camacho
  1489. * Changes to return discontinuous seq-aligns for each query-subject pair
  1490. *
  1491. * Revision 1.20  2003/10/30 21:40:57  dondosha
  1492. * Removed unneeded extra argument from BLAST_Results2CSeqAlign
  1493. *
  1494. * Revision 1.19  2003/10/28 20:53:59  camacho
  1495. * Temporarily use exception for unimplemented function
  1496. *
  1497. * Revision 1.18  2003/10/15 16:59:42  coulouri
  1498. * type correctness fixes
  1499. *
  1500. * Revision 1.17  2003/09/09 15:18:02  camacho
  1501. * Fix includes
  1502. *
  1503. * Revision 1.16  2003/08/19 20:27:51  dondosha
  1504. * Rewrote the results-to-seqalign conversion slightly
  1505. *
  1506. * Revision 1.15  2003/08/19 13:46:13  dicuccio
  1507. * Added 'USING_SCOPE(objects)' to .cpp files for ease of reading implementation.
  1508. *
  1509. * Revision 1.14  2003/08/18 22:17:36  camacho
  1510. * Renaming of SSeqLoc members
  1511. *
  1512. * Revision 1.13  2003/08/18 20:58:57  camacho
  1513. * Added blast namespace, removed *__.hpp includes
  1514. *
  1515. * Revision 1.12  2003/08/18 17:07:41  camacho
  1516. * Introduce new SSeqLoc structure (replaces pair<CSeq_loc, CScope>).
  1517. * Change in function to read seqlocs from files.
  1518. *
  1519. * Revision 1.11  2003/08/15 15:56:32  dondosha
  1520. * Corrections in implementation of results-to-seqalign function
  1521. *
  1522. * Revision 1.10  2003/08/12 19:19:34  dondosha
  1523. * Use TSeqLocVector type
  1524. *
  1525. * Revision 1.9  2003/08/11 14:00:41  dicuccio
  1526. * Indenting changes.  Fixed use of C++ namespaces (USING_SCOPE(objects) inside of
  1527. * BEGIN_NCBI_SCOPE block)
  1528. *
  1529. * Revision 1.8  2003/08/08 19:43:07  dicuccio
  1530. * Compilation fixes: #include file rearrangement; fixed use of 'list' and
  1531. * 'vector' as variable names; fixed missing ostrea<< for __int64
  1532. *
  1533. * Revision 1.7  2003/08/01 17:40:56  dondosha
  1534. * Use renamed functions and structures from local blastkar.h
  1535. *
  1536. * Revision 1.6  2003/07/31 19:45:33  camacho
  1537. * Eliminate Ptr notation
  1538. *
  1539. * Revision 1.5  2003/07/25 22:12:46  camacho
  1540. * Use BLAST Sequence Source to retrieve sequence identifier
  1541. *
  1542. * Revision 1.4  2003/07/25 13:55:58  camacho
  1543. * Removed unnecessary #includes
  1544. *
  1545. * Revision 1.3  2003/07/23 21:28:23  camacho
  1546. * Use new local gapinfo structures
  1547. *
  1548. * Revision 1.2  2003/07/14 21:40:22  camacho
  1549. * Pacify compiler warnings
  1550. *
  1551. * Revision 1.1  2003/07/10 18:34:19  camacho
  1552. * Initial revision
  1553. *
  1554. *
  1555. * ===========================================================================
  1556. */