alngraphic.cpp

来自「ncbi源码」· C++ 代码 · 共 696 行 · 第 1/2 页

/*
 * ===========================================================================
 * PRODUCTION $Log: alngraphic.cpp,v $
 * PRODUCTION Revision 1000.0  2004/06/01 19:48:35  gouriano
 * PRODUCTION PRODUCTION: IMPORTED [GCC34_MSVC7] Dev-tree R1.2
 * PRODUCTION
 * ===========================================================================
 */

/*  $Id: alngraphic.cpp,v 1000.0 2004/06/01 19:48:35 gouriano Exp $
 * ===========================================================================
 *
 *                            PUBLIC DOMAIN NOTICE
 *               National Center for Biotechnology Information
 *
 *  This software/database is a "United States Government Work" under the
 *  terms of the United States Copyright Act.  It was written as part of
 *  the author's official duties as a United States Government employee and
 *  thus cannot be copyrighted.  This software/database is freely available
 *  to the public for use. The National Library of Medicine and the U.S.
 *  Government have not placed any restriction on its use or reproduction.
 *
 *  Although all reasonable efforts have been taken to ensure the accuracy
 *  and reliability of the software and data, the NLM and the U.S.
 *  Government do not and cannot warrant the performance or results that
 *  may be obtained by using this software or data. The NLM and the U.S.
 *  Government disclaim all warranties, express or implied, including
 *  warranties of performance, merchantability or fitness for any particular
 *  purpose.
 *
 *  Please cite the author in any work or product based on this material.
 *
 * ===========================================================================
 *
 * Author:  Jian Ye
 *
 * File Description:
 *   Alignment graphic overview (using HTML table) 
 *
 */

#include <ncbi_pch.hpp>
#include <objtools/align/alngraphic.hpp>
#include <util/range.hpp>
#include <serial/iterator.hpp>
#include <objects/seq/Bioseq.hpp>
#include <objects/seqloc/Seq_loc.hpp>
#include <objects/seqloc/Seq_id.hpp>
#include <objects/seqalign/Seq_align.hpp>
#include <objects/seqalign/Seq_align_set.hpp>
#include <objects/seqalign/Score.hpp>
#include <objects/seqalign/Dense_seg.hpp>
#include <objects/seqalign/Std_seg.hpp>
#include <objects/seqalign/Dense_diag.hpp>
#include <objects/seqalign/seqalign_exception.hpp>
#include <objects/general/Object_id.hpp>
#include <objmgr/util/sequence.hpp>
#include <html/html.hpp>

BEGIN_NCBI_SCOPE
BEGIN_SCOPE (objects)
USING_SCOPE (sequence);

static string kDigitGif[] = {"0.gif", "1.gif", "2.gif", "3.gif", "4.gif",
                             "5.gif", "6.gif", "7.gif", "8.gif", "9.gif"};

static const TSeqPos kScoreMargin = 50;
static const TSeqPos kScoreLength = 500;
static const TSeqPos kScoreHeight = 40;

static const TSeqPos kBlankBarHeight = 4; 

static const TSeqPos kMasterHeight = 10;
static const TSeqPos kMasterPixel = kScoreLength;
static const TSeqPos kMasterBarLength = 550;

static const TSeqPos kScaleMarginAdj = 1;
static const TSeqPos kScaleWidth = 2;
static const TSeqPos kScaleHeight = 10;

static const TSeqPos kDigitWidth = 10;
static const TSeqPos kDigitHeight = 13;

static const TSeqPos kGapHeight = 1;

static const TSeqPos kNumMark = 6;

static const int kDeflineLength = 55;

/*two ranges are considered overlapping even if they do not overlap but
  they differ only by this number*/
static const int kOverlapDiff = 5;

//gif images
static const string kGifWhite = "white.gif";
static const string kGifMaster ="query_no_scale.gif";
static const string kGifScore = "score.gif";
static const string kGifGrey = "grey.gif";
static const string kGifScale = "scale.gif";

static string s_GetGif(int bits){
    string gif = NcbiEmptyString;
    if(bits < 40){
        gif = "black.gif";
    } else if (bits < 50) {
        gif = "blue.gif";
    } else if (bits < 80) {
        gif = "green.gif";
    } else if (bits < 200) {
        gif = "purple.gif";
    } else  {
        gif = "red.gif";
    } 
    //for alignment with no bits info
    if(bits == 0){
        gif = "red.gif";
    }
    return gif;
}

//Round the number.  i.e., 290 => 250, 735 => 700. 
static int s_GetRoundNumber (int number){
    int round = 0;
   
    if (number > 10) {
        string num_str = NStr::IntToString(number);
        int trail_num = NStr::StringToInt(num_str.substr(1, num_str.size() - 1));
        round = number - trail_num;
        string mid_num_str;
        int mid_num = 0;
        if(num_str.size() > 2) {
            string trail_zero (num_str.size() - 2, '0');
            mid_num_str = "5" +  trail_zero;
            mid_num = NStr::StringToInt(mid_num_str);
        }
       
        if (number >= round + mid_num) {
            round += mid_num;
        } 
    } else {
        round = number;
    }
    return round;
}

CRange<TSeqPos>* CAlnGraphic::x_GetEffectiveRange(TAlnInfoList& alninfo_list){
    CRange<TSeqPos>* range = NULL;
    if(!alninfo_list.empty()){
        range = new CRange<TSeqPos>(alninfo_list.front()->range->GetFrom(), 
                                    alninfo_list.back()->range->GetTo());
    }
    return range;
}

void CAlnGraphic::x_MergeDifferentSeq(double pixel_factor){
  
    TAlnInfoListList::iterator iter_temp;
    list<CRange<TSeqPos>* > effective_range_list;
    CRange<TSeqPos> *effective_range = NULL, *temp_range = NULL;
   
    NON_CONST_ITERATE(TAlnInfoListList, iter, m_AlninfoListList) {
        iter_temp = iter; //current list
        iter_temp ++;  //the next list

        if(!(*iter)->empty() && iter_temp != m_AlninfoListList.end()){
            /*effective ranges means the range covering all ranges for
              alignment from the same seq*/
            temp_range = x_GetEffectiveRange(**iter);
            if(temp_range){
                effective_range_list.push_back(temp_range);
            }
  
            /*compare the range in current list to range in all other lists.
              If the latter does not overlap with the ranges in the current
              list, move it to the current one*/
            while(iter_temp != m_AlninfoListList.end()){
                /* Get effective range. Note that, except the current list,
                   all other list contains only one effective range. The
                   current list may have more than one because the range 
                   from other list may be moved to it*/
                effective_range = x_GetEffectiveRange(**iter_temp);
                bool overlap = false;
                if(effective_range){
                    NON_CONST_ITERATE(list<CRange<TSeqPos>*>,
                                      iter2, 
                                      effective_range_list){
                        CRange<TSeqPos>* temp_range2 = *iter2;
                        int min_to = min(effective_range->GetTo(), 
                                         temp_range2->GetTo());
                        int max_from =  max(effective_range->GetFrom(), 
                                            temp_range2->GetFrom());
                        //iter2 = ranges in current list
                        if((int)((max_from - min_to)*pixel_factor) <= kOverlapDiff){
                            //Not overlap but has only overlap_diff is still
                            //considered overlap as it's difficult to see the
                            //difference of only 1 pixel 

                            overlap = true;
                            break;  //overlaps, do nothing
                        } 
                    }
                    //if no overlap, move this list to the current list
                    if(!overlap){
                        (*iter)->merge(**iter_temp);
                        effective_range_list.push_back(effective_range);
                    }
                    if(overlap){
                        delete effective_range;
                    }
                }
                iter_temp ++;
            }
            ITERATE(list<CRange<TSeqPos>* >, iter4, effective_range_list){
                delete *iter4;
            }
            effective_range_list.clear();
        }
        (*iter)->sort(FromRangeAscendingSort);
    }
}
void CAlnGraphic::x_MergeSameSeq(TAlnInfoList& alninfo_list){  
    TAlnInfoList::iterator prev_iter; 
    int i = 0;
 
    NON_CONST_ITERATE (TAlnInfoList, iter, alninfo_list){
        if(i > 0 && (*prev_iter)->range->IntersectingWith(*((*iter)->range))){
            //merge range
            (*iter)->range->Set(min((*prev_iter)->range->GetFrom(), 
                                    (*iter)->range->GetFrom()),
                                max((*prev_iter)->range->GetTo(), 
                                    (*iter)->range->GetTo()));
            //take the maximal score and evalue 
            if((*iter)->bits < (*prev_iter)->bits){
                (*iter)->bits = (*prev_iter)->bits;
                (*iter)->info = (*prev_iter)->info;
            }
            delete (*prev_iter)->range;
            delete *prev_iter;
            alninfo_list.erase(prev_iter);     
        }
        i ++;
        prev_iter = iter;
    }    
}

template<class container> 
static bool s_GetBlastScore(const container& scoreList,  
                                                      int& score, 
                                                      double& bits, 
                                                      double& evalue){
    bool hasScore = false;
    ITERATE (typename container, iter, scoreList) {
        const CObject_id& id=(*iter)->GetId();
        if (id.IsStr()) {
            hasScore = true;
            if (id.GetStr()=="score"){
                score = (*iter)->GetValue().GetInt();

            } else if (id.GetStr()=="bit_score"){
                bits = (*iter)->GetValue().GetReal();

            } else if (id.GetStr()=="e_value" || id.GetStr()=="sum_e") {
                evalue = (*iter)->GetValue().GetReal();
            } 
        }
    }
    return hasScore;
}

static void s_GetAlnScores(const CSeq_align& aln, 
                           int& score, 
                           double& bits, 
                           double& evalue){
    bool hasScore = false;
    //look for scores at seqalign level first
    hasScore = s_GetBlastScore(aln.GetScore(),  score, bits, evalue);
    
    //look at the seg level
    if(!hasScore){
        const CSeq_align::TSegs& seg = aln.GetSegs();
        if(seg.Which() == CSeq_align::C_Segs::e_Std){
            s_GetBlastScore(seg.GetStd().front()->GetScores(), 
                            score, bits, evalue);
        } else if (seg.Which() == CSeq_align::C_Segs::e_Dendiag){
            s_GetBlastScore(seg.GetDendiag().front()->GetScores(), 
                            score, bits, evalue);
        }  else if (seg.Which() == CSeq_align::C_Segs::e_Denseg){
            s_GetBlastScore(seg.GetDenseg().GetScores(),  score, bits, evalue);
        }
    }
}


void CAlnGraphic::x_GetAlnInfo(const CSeq_align& aln, const CSeq_id& id,
                               SAlignInfo* aln_info){
    string info;
    int score = 0;
    double evalue = 0;
    double bits = 0;
    string bit_str, evalue_str;
    string title;
  
    const CBioseq_Handle& handle = m_Scope->GetBioseqHandle(id);
    if(handle){
        const CBioseq::TId& ids = handle.GetBioseqCore()->GetId();
        CRef<CSeq_id> wid = FindBestChoice(ids, CSeq_id::WorstRank);
        aln_info->id = wid;
        aln_info->gi =  FindGi(ids);
        wid->GetLabel(&info, CSeq_id::eContent, 0);
        title = GetTitle(handle);

        if (title.size() > kDeflineLength){
            title = title.substr(0, kDeflineLength) + "..";
        }
        info += " " + title;   
    } else {
        aln_info->gi = 0;
        aln_info->id = &id;
        aln_info->id->GetLabel(&info, CSeq_id::eContent, 0);
    }
    s_GetAlnScores(aln, score, bits, evalue);
    NStr::DoubleToString(bit_str, (int)bits);
    NStr::DoubleToString(evalue_str, evalue);
    
    CNcbiOstrstream ostream; 
    ostream << setprecision (2) << evalue;
    string formatted_evalue = CNcbiOstrstreamToString(ostream);
    info += " S=" + bit_str + " E=" + formatted_evalue;
    aln_info->info = info;
    aln_info->bits = bits;
   
}

CAlnGraphic::CAlnGraphic(const CSeq_align_set& seqalign, CScope& scope)
    :m_AlnSet(&seqalign), m_Scope(&scope)  {
    m_NumAlignToShow = 1200;
    m_View = eCompactView;
    m_BarHeight = e_Height4;
    m_ImagePath = "./";
    m_MouseOverFormName = "document.forms[0]";
    m_NumLine = 55;
}

CAlnGraphic::~CAlnGraphic(){
    //deallocate memory
    ITERATE(TAlnInfoListList, iter, m_AlninfoListList) {
        ITERATE(TAlnInfoList, iter2, **iter){
            delete (*iter2)->range;